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INTRODUCTION 

Spectral  analysis  of  stationary  random  processes  via  linear  predictive, 
maximum  entropy,  and  autoregressive  techniques  has  attracted  much  attention 
lately,  especially  for  short  data  segments;  see,  for  example,  the  biblio- 
graphies listed  in  references!,  2,  and  3.  For  a univariate  process,  it 
appears  that  the  Burg  algorithm  (Ref.  4),  which  guarantees  a stable  correla- 
tion recursion,  is  as  good  as  any  of  the  currently  available  techniques  of 
similar  nature  that  employ  an  all-pole  model  of  the  available  process  (Ref.  3). 

Accordingly,  it  is  desirable  to  develop  a spectral  analysis  technique 
for  the  multivariate  case  in  such  a way  that:  we  employ  a physically  mean- 
ingful error  minimization  for  the  determination  of  the  filter  coefficients; 
the  technique  yields  a stable  correlation  recursion;  and  it  reduces  to  Burg's 
algorithm  for  the  univariate  case.  It  will  be  shown  in  the  following  that 
we  have  accompl i shed  these  goals,  with  the  exception  that  we  have  not  proved 
(or  disproved)  the  stability  requirement.  A FORTRAN  program  for  this  spectral 
analysis  technique  was  published  in  Ref.  5,  along  with  an  example  of  its 
application.  Virtually  simultaneously,  the  same  technique  was  investigated 
independently  and  published  in  Ref.  6.  In  this  report,  we  will  document  the 
derivations  and  equations  that  lead  to  the  program  presented  in  Ref.  5,  and 
indicate  an  extension  of  that  result. 

Our  approach  in  this  report  will  be  to  investigate,  in  some  detail,  first 
the  case  where  the  correlation  of  the  multivariate  process  under  consideration 
is  known  for  a limited  range  of  argument  values,  and  to  extract  all  the 
relevant  important  properties  of  the  solution  so  that  they  may  be  forced  to 
be  satisfied  later  when  we  treat  the  unknown  correlation  case.  This  property- 
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extraction  procedure  will  be  found  to:  furnish  guides  to  the  analysis  of 
the  unknown  correlation  case;  allow  us  to  cut  down  on  computer  execution 
time  and  storage  by  employing  the  properties;  and  make  us  aware  of  some  of 
the  shortcomings  of  the  unknown  (versus  known)  correlation  cases.  This 
procedure  should  also  be  helpful  to  those  who  are  not  thoroughly  familiar 
with  spectral  analysis  of  multivariate  processes  and  their  properties. 

Throughout  this  report,  we  assume  we  are  dealing  with  equi spaced 
samples  of  a stationary  zero-mean  complex  random  process  X(t)  of  dimension- 
ality M;  that  is,  sample 

xW»x«[<-<’]T  <’) 

is  an  M x 1 column  matrix,  where  A is  the  common  sampling  interval  for  all 
the  component  processes  of  X(t).  It  is  not  assumed  that  X(t)  is  Gaussian. 

In  section  2,  we  will  assume  that  the  correlation  matrix  of  process 
(Xn),  namely  the  M x M matrix* 

K, ■ xjC  * (2) 

is  known  exactly  for  a limited  range  of  values  of  k,  and  will  show  how 
an  approximation  for  tne  spectrum  of  process  {Xn}  can  be  obtained.  In 
section  3,  the  input  correlation  matrix  R*  will  be  unknown,  and  all  that 
is  available  is  a finite  set  of  N data  samples,  X,,  x2,  ...,  XN,  from 
which  an  estimate  of  the  spectrum  of  process  { Xn}  is  desired.  The  end 
result  will  be  a FORTRAN  program  for  multivariate  spectral  analysis. 

*The  case  of  complex  samples  is  treated  so  that  we  can  handle  complex 
envelope  or  complex  demodulated  processes.  Specialization  to  real  processes 
is  immediate,  and  (2)  becomes  R„=R^.  An  overbar  indicates  an  ensemble 
average,  superscript  T denotes  a transpose,  and  superscript  H denotes  a 
conjugate  transpose.  Matrices  are  indicated  by  capital  letters. 
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2.  KNOWN  CORRELATION 


If  the  correlation  in  (2)  is  known  for  all  k,  the  standard  (double- 
sided) definition  of  the  spectrum  of  process  (Xn)  is 

6c  (f)  * A J^xp(-i'brfkA)  R*  , lf)<jk  • ^ 

The  complex  M x M matrix  6(f)  is  Hermitian  and  non-negative  definite  for 
any  value  of  frequency  f (see  appendix  A),  but  need  not  be  even  in  frequency 
f.  When  Rj,  is  not  known  for  all  k,  but  only  for  a range  an  approxi- 

mation to  (3)  must  bo  accepted;  this  problem  will  be  pursued  below. 


2.1  DERIVATION  OF  EQUATIONS 


Suppose  M-dimensional  samples  X.  , ...,  X.  are  available,  and  we 

K“P  K- » 

attempt  a one-step  linear  prediction  of  X^  according  to  the  p-th  order 


operation 


where  complex  coefficient  matrix  Ap  is  M x M,  n = 1 , 2,  p.  The 
instantaneous  error  at  time  kA  is  defined  as 

x*  x-x,  Iax.  - A’*~T-  (bi 

The  linear  operators  in  (4)  and  (5)  constitute  stable  linear  filters  regard- 
less of  the  choice  of  coefficients;  the  filter  of  (4)  is  called  the  predictive 
filter,  that  of  (5)  is  called  the  predictive  error  filter.  Notice  that  we  are 
not  assuming  that  process  { Xn } actually  satisfies  an  autoregressive  relation; 
rather  we  are  simply  attempting  to  linearly  predict  { Xp } on  the  basis  of  the 
most  recent  p past  values. 
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The  minimum  value  of  the  scale*,  error 


X“X  • ir  \\H  , 

p 


(6) 


by  choice  of  coefficients  {Anl 2 , is  given  (in  appendix  B)  by  the  solution 
of  the  linear  matrix  equations 

(7) 

where  the  explicit  dependence  on  the  order  p is  indicated.  Knowledge  of 
for  Ikkp  is  required  in  (7). 

Before  we  discuss  the  solution  of  (7)  for  , we  consider  one-step 

linear  "backward  prediction"  of  process  (Xn).  Suppose  samples  Xk,  Xk_i, 
Xk_p+i  are  available,  and  we  attempt  a one-step  linear  prediction  of  Xk_p 
according  to 

V ' (8) 

The  instantaneous  error  is  defined  as 


, V-r- 


(9) 


(10) 


The  minimum  value  of  the  scalar  error 

by  choice  of  coefficients  (Bn)? , may  be  shown  (in  a manner  similar  to  that 
of  appendix  B)  to  be  given  by  the  solution  of  the  linear  matrix  equations 

'"Kn  > J - f-  (11) 

For  the  optimum  coefficients  in  (7)  and  (11),  we  find  (see  appendix  B) 
that  the  optimum  error  matrices  take  the  form 
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oft  xx”  = K-  £a!il  = Op , U*  R. , 

.rt^r  = K-iBTR.'  % , Y-=  K. 


02) 


In  general,  these  two  matrices,  their  diagonal  elements,  and  their  traces  are 
unequal  (as  the  simple  example  of  p=l  will  show).  However,  their  determinants 
are  equal,  as  will  be  shown  in  subsection  2.2. 

The  solutions  of  (7)  and  (11)  can  be  accomplished  simultaneously  in  a 
recursive  fashion  (Ref.  7).  Define 

A?1; -I, 

pp..-£bTV,  , Bf-I. 


03) 


Then 

and 


—i 


< * cr  v;; , bjt-  ^ 


Jr0 


AT  - A.  ty. 


jF-bT-eCA^ 


(14) 

05) 


| < M < }>-»  (p  - 2)  • 

r np-" 

These  relations  will  be  simplified  somewhat  in  subsection  2.2.  For  M=l, 
a univariate  process,  (7)  and  (11)  immediately  yield 

* - 06) 


a!’-bT  £,r  M-l, 


where  we  have  used  (2)  in  the  form  R = R for  a univariate  process.  No 

k -k 


such  simple  relation  as  (16)  holds  for  M > 2. 


We  will  now  derive  a chain  interpretation  of  the  above  results  that  will 


prove  very  useful  later  when  we  have  to  deal  with  the  unknown  correlation  case 

P . P 


\ , define  the  p-th 

i i;nJ,  K 


For  the  optimum  filter  coefficients  ( ' and 
order  forward  and  backward  residuals  (see  (5)  and  (9))  as  the  outputs  of  the 
forward  and  backward  predictive  error  filters: 


TR  5501 


^’,-iA!xK.„-X-AfXw--ArX.f  . 

* W'O 


Then  using  (15),  we  can  expre^ 


And  similarly 


X^V|aX-<x,, 

= K-KAr-ArcK-^x, 


. vrv  a!’  2tf*X  = Y*~Ar>2<M 

* j=0  J ‘k  nf  • 


2f- V £*X.*bTx. 


(17) 


(18) 


(19) 


\/^p)  “7  (p) 

Thus  p-th  order  residuals  YK  and£^  are  related  to  the  ( p-1 ) th  order  residuals 
simply  through  the  coefficients  A'jf  and  5‘f.  A block  diagram  of  the  relation- 
ships in  (18)  and  (19)  is  given  in  figure  1,  where  2^1  denotes  an  M x M matrix 


filter  of  unit  delay. 
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-(P-1) 


(P-D| 


- 1 T \ (O)  (O) 

2 1 )Yl  =Z„  = 

v y k k k 


Figure  1.  Chain  Representation  of  Residuals 
Thus  matrix  operators  ^ and  can  be  interpreted  as  those  coefficients 


which  minimize 


r wnP”' 

respectively,  at  the  output  of  the  p-th  stage  in  figure  1,  where^  J and 
«rr  are  determined  by  minimizations  at  lower  order  stages.  A**  and  8^ 
are  called  the  partial  correlation  coefficients.  Stated  alternatively,  stage 
by  stage  minimizations  of  (20),  via  choices  of  partial  correlation  coefficients 
and  respectively,  results  in  the  same  overall  filter  as  if  the  powers 


- i/tx 


were  minimized  by  the  choices  of  (An)?  and  {Bn }? , respectively,  each  in  one 
simultaneous  optimization.  This  will  furnish  an  important  reference  point 
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for  the  unknown  correlation  case  in  section  3. 

If  we  let  the  transfer  functions  (z- transforms)  to  the  outputs  of  the 
p-th  stage  in  figure  I be  denoted  by^(z)  and  it  immediately 

follows,  from  figure  1 or  equations  (18)  and  (13),  that 


ip> 

► 

A 

¥ 


(*)  - 
w* 


> 

*-*rv  Bjr^rw, 


= x . 


(22) 


In  closed  form,  these  predictive  error  filter  transfer  functions  are  express- 
ible  as  (see  (17))  = - £ g ’ g , j . 

* Mc«  " W-l  '*»  ) 


u w*l  J • 


(23) 


2.2  PROPERTIES  AND  INTERPRETATIONS 


Suppose  that  process  { Xn > were  scaled  according  to 

X*  - T>  X. 

where  M x M matrix  D is  arbitrary,  but  invertible.  Then  the  correlation  of 
the  scaled  process  is 


(24) 


\ * X )£  * xJ^d”  = dk,.xm. 

Now  from  (7),  since  the  solutions  ul)  and  m must  satisfy 

^ ^ r ft*  Is  > 

*■ ) ) ' 


(25) 


tOk\ 

\ l 
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? 'A 


\ i 


t 

i 


.V  ? 

t 1 

? i 


! I 


I > 


•?  i 

f * 


V 


^AfK,.„=t, , i*”sp> 


(27) 


Vt  = I 


respectively,  the  solutions  are  related  by  a similarity  transformation: 

AM>A!V,  i.Hip.  (28) 

This  is  called  the  scaling  property.  A similar  property  holds  for  the  back- 
ward coefficients  Wl. 


An  immediate  by-product  of  the  scaling  property  is  that  A^  and  A have 


the  same  eigenvalues: 


U{hP  XT)  = Xl)  = cki (a!’- X r) . 


(29) 


Similarly,  and  B?  have  the  same  eigenvalues,  regardless  of  scaling  matrix 


The  remainder  of  this  subsection  will  deal  .vith  the  quantities  Up  and 
Vp  defined  in  (12),  and  Cp  and  Dp  defined  in  (13).  The  quantity  Up  can  be 
interpreted  physically  as  the  correlation  matrix  of  the  p-th  order  forward 
residual;  see  (12),  (5),  and  (17).  Similarly,  Vp  is  the  correlation  matrix 
of  the  p-th  order  backward  residual;  see  (12),  (9),  and  (17).  That  is, 


Ur = V XT"  , \ 


(30) 


Thus  U and  V are  Hermitian: 
P P 


(31) 


and  Up  and  V are  non-negative  definite: 


= y H'ATT- 0 


(32) 


for  any  M x 1 matrix  "f  . In  appendix  C,  it  is  shown  that  simple  recur- 


sions hold  for  Up  and  Vp: 


5 


'i  ’ 

1; 

i ' 


Tv 


£ 
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Ji 
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0,.  (I-S>')uri  , Ll-t, 

V(i-tF <)  V,., , V ■ K. . 

It  immediately  follows  from  (33)  that  (see  appendix  C) 

<W  Uf-  4etVf  ; p&o. 

This  property  was  proved  in  Ref.  8,  page  240. 

The  quantity  Cp  defined  in  (13)  can  be  interpreted  as  the  cross-corre- 
lation matrix  between  the  p-th  order  forward  and  backward  residuals  at  one 
unit  of  delay:  rfrZfiX  A.  L.  TT  h~ 

* Z«-'  = ^ * Ak-n  , 

where  we  have  used  (17),  (2),  (7),  and  (13),  Similarly 

where  we  have  used  (17),  (2),  (11),  and  (13).  It  immediately  follows  from 
(35)  and  (36)  that 

v?  - S . , 

Thus  it  is  not  necessary  to  do  the  additional  calculation  of  Dp  in  the 
solution  given  in  (13). 

Another  interpretation  of  Cp  is  available  as  follows: 

tr  .^~£e>uZ-  - i«*.-  ■(»,',=-  p|( 

t S' 

where  we  have  employed  (17),  (2),  (12),  (7),  and  (13)  in  order.  Thus  the 
p-th  order  forward  residual  is  uncorrelated  with  the  p most  recent  past 


• - o 
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values  of  the  input,  and  the  crosscorrelation  at  p + 1 units  of  delay  is 
just  Cp.  Similarly,  the  backward  residual  satisfies 


Vf , *1=0 


s = ] 0>  »1W,1PL(39) 

Dpy  W " p-H 


Yet  another  interpretation  of  Cp  and  Dp  will  be  given  in  subsection  2.3. 
As  the  order  p in  the  linear  prediction  (4)  increa:s?  (38)  yields 


V H 


u /m=oi 


t,  U Y*\ 


Q5  j>“*  ^ 


(40) 


Therefore  the  autocorrelation  matrix  of  the  forward  residual  becomes 


*s  °°-  (41) 


That  is,  p-th  order  residual  Yjf  tends  to  white  noise  with  a correlation 
matrix  at  zero  time  delay  of  value  Uro,  which  is  not  necessarily  diagonal 

The  Herrnitian  property  in  (37)  allows  us  to  combine  (14)  into  the 
equation 

[V 


a;\,  - ^ , 


(42) 
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where  we  utilized  (31).  This  constraint  on  the  partial  correlation  coeffi- 
cients will  be  of  paramount  importance  in  the  unknown  correlation  case.  It 
immediately  follows  from  (42)  and  (34)  r.hat 


No  such  simple  relation  holds  between  det  and  det  B^for  n<p,  except 
for  M = 1 , a univariate  process. 

2.3  EXTRAPOLATION  OF  CORRELATION  VALUES 

In  subsection  2.1,  we  minimized  the  error  in  prediction  (4)  and  found 
that  for  a p-th  order  prediction,  knowledge  of  Rk  for  IKS* p was  re  ired; 

see  (7).  Now  suppose  that  this  is  all  the  knowledge  available  about  iRk) ; 
that  is,  suppose  Rk  i s unknown  for  jk|^p.  What  can  be  done  about  approximating 
these  unknown  values? 

One  approach  is  as  follows:  we  assume  that  the  p-th  order  residual 
process  in  (17)  is  white  (i.e.,  uncorrelated  for  all  non-zero  delays),  and 

that  /^=o(otherwise  we  could  reduce  the  value  of  p).  That  is,  we  assume  we 
can  do  nothing  more  in  prediction  by  choosing  more  terms  in  the  sum  (4), 
which  is  tantamount  to  assuming  maximum  uncertainty  (entropy)  about  the 
residual  process  jY„r}-  This  is  a very  extensive  assumption;  wo  now  investi- 
gate its  ramifications. 

We  know  from  (38)  that 


12 


must  satisfy 


^ = 0 | :£  m 5 t>. 


( 


Additionally,  employing  (17),  the  autocorrelation  matrix  of  the  p-th  order 
residual  is 


( 


Now  for  j = 1,  the  white  noise  assumption  on  process  {Nff’}yields , via  (46) 
and  (45), 

«--*F.w--E?xr;u  £.0.  . 


'*»+!  ' '»  f>+l  ' *P  ) ’ff* 

And  for  j = 2,  the  white  assumption  (in  conjunction  with  (47))  yields 


tf-o. 


W“-U 


Continuing  in  this  way,  the  white  assumption  is  tantamount  to  assuming  that 

- 0 "for  j>-H  S yv).  ^ 


Returning  to  expression  (44),  this  means  that  we  are  assuming  that 

* 0 phs  m>  (5 

that  is, 

^ pf/  s K»l.  (5 


Using  more  explicit  notation,  and  denoting  these  assumed  values  of  correlati 
as  forward  extrapolations  , we  have 
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p + l * W, 


(52) 


where  "starting  values" 

R»  * } 0 5 * < j>.  (53) 

Equation  (52)  is  called  the  correlation  recursion  equation.  It  is  interesting 
to  note  that  the  form  of  the  correlation  recursion  (52)  is  identical  to  the 
form  (4)  for  the  individually  predicted  waveform  values. 


The  correlation  values  in  (52)  are  called  the  maximum  entropy  correlation 
extrapolations.  The  recursion  is  stable  if  and  only  if  (see  (23)) 

*-/£)•  (54) 

' If®  ) ''  ' 

possesses  all  its  zeros  within  the  unit  circle  in  the  complex  z-plane; 
this  property  will  be  treated  in  subsection  2.4. 


A similar  procedure  for  backward  correlation  extrapolation,  assuming 
that  residual  process  [£**]  is  white,  yields 


fc'l  V, 


r 


•H  — W 


(55) 


where 


0 ^ wt  < 


T 


(5b) 
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Backward  recursion  (55)  is  identical  in  form  to  the  backward  prediction  (8). 
The  recursion  (55)  is  stable  if  and  only  if  (see  (23)). 

det(T-£2-"Bt)=  (57) 

possesses  all  its  zeros  within  the  unit  circle. 


As  a special  case  oi  (52)  and  (53),  the  one-step  forward  extrapolated 
correlation  based  on  a p-th  order  prediction  is 


r «»-i  r*  w=  \ J 


But  from  (13),  we  now  can  see  that 


p+v-» 


v» * O 


(58) 


(59) 


That  is,  Cp  is  the  difference  between  the  true  correlation  value  Rp+]  and  the 
one-step  forward  extrapolated  correlation  Rp+1  based  upon  knowledge  of  { R^ } P 


A similar  procedure  shows  that 

D„  - 1? 


f ■ v.  - r;:  . 


p. 


(60) 


That  is.  Dp  is  the  difference  between  the  true  correlation  value  R_p_-|  and 


the  one-step  backward  extrapolated  correlation  based  upon  knowledge 


°f  ^k^  • 


When  (59)  and  (60)  rp  combined  with  the  Hermitian  property  in  (37),  we 
see  that 


(61) 


JK 

•J 


$ 
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This  is  a special  case  of  the  mere  general  property  (demonstrated  in  appendix 
D)  that 

’ft-T  * Klf  , f + l s ">  (62) 


that  is,  the  backward  and  forward  extrapolated  correlation  matrices  are 
Hermitians  of  each  other.  This  is  a desirable  property  of  the  extrapolations 
and  is  consistent  with  the  same  property,  (2),  which  holds  for  the  known 

p 

correlation  values, 


It  was  noted  in  (54)  and  (57)  that  the  zeros  of  detft^(z)  and  det'ftjfU"1 ) 
must  be  within  the  unit  circle  in  order  that  recursions  (52)  and  (55),  respect- 
ively, be  stable.  It  is  shown  in  appendix  E that 


det  (i  • - ± «- M?’ ) - to  (i  ■ - i i*  bVh  ) 

Mr  J ' ' >•-)  / ' 


(63) 


That  we  need  consider  only  the  zeros  of  one  of  thtse  quantities;  the  location 
of  these  zeros  is  considered  below. 


It  is  also  shown  in  appendix  E that 

tr  a*  . MT)* 

and 

aa  a*/  = (aa  . 

2.4  SPECTRAL  APPROXIMATION 


(64) 


(65) 


| Equations  (52)  and  (53)  define  the  forward  extrapolated  correlations  for 

j Ail  m>0.  We  extend  these  to  negati  /e  m via 

I 


r1:’  * k 


S 0. 


(66) 
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which  is  consistent  with  (2).  We  will  now  use  the  Fourier  transform  of  this 
infinite  sequence,  as  in  (3),  as  an  approximation  to  the  spectrum  of  process 
{Xn}.  In  appendix  F,  it  is  shown  that  the  approximate  spectrum  is  given  by 


where 


<?(A  - - W Up  ft"  , If)  < i * 

hIw  * - i:  ••2rf-*)4T 


is  the  forward  predictive  error  filter  transfer  function.  Since  Up  is  non- 
negative  definite  by  (32),  spectral  approximation  Gvp;(f)  is  nonnegative 
definite  for  any  f;  it  is  also  obviously  Hermitian  by  (31).  Thus  the  desir- 
able properties  of  appendix  A are  achieved  by  approximation  (67).  In  order 
to  evaluate  (67),  one  M x M matrix  inverse  (of  H^(f))  is  needed  at  each 
value  of  f of  interest. 


A similar  pocedure  applied  to  the  backward  correlation  recursion  of 
(55)  and  (56)  yields  the  spectral  approximation 


where 


G?(f)  -4  Hftf  VP  Hf(f)"  , lfl«i 

^ e<p(-’2-rrf«a)3^, 


is  the  backward  predictive  error  filter  transfer  function.  Since  the  extra- 
polated correlations  via  (52)  or  (55)  are  equal,  as  shown  in  subsection  2.3, 
the  same  notation  G^’)(f ) , is  used  for  both  (67)  and  (69);  however,  we  have 
two  different  factorizations  for  the  unique  spectral  approximation  G^P^(f). 

In  appendix  F,  it  is  also  shown  that  the  zeros  of  det^(r)  (see  (22) 
and  (23))  all  lie  inside  the  unit  circle  in  the  complex  z-plane.  Additionally, 
the  poles  of  ' all  lie  inside  the  unit  circle,  and  the  zeros  offtMj.)  1 
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? i 

? 1 


all  lie  at  z = 0.  Thus  the  recursion  (52)  is  stable.  This  point  is 
discussed  in  Ref.  7,  p.  132. 


2.5  EXAMPLE 


A simple  example  for  M = 2 will  be  considered.  Let  the  process  be 
generated  according  to 


where 


Q - 


.25  -.15- 

65  . Sir 


and  white  noise  satisfies 


W * LJ- 

Then  it  may  be  shown  that 

K, * h r . *-0, 


with  solution 


)35 
j , 

4.  K2 

Tf,  = 

nr ns 

-12. 

L 4.K2 

21-C4  3_ 

/ 1 

no  12 

If.  064 

and  (11), 

we  find 

J 

k- 


.85 

-.is7 

4°- 

.55*130 

.7 5279 

.65 

.55^ 

> 

i 

(>44©0 

. ?4070 

and  A, -A,  o,  2 S n < |>  . We  observe 

and  The  determinants  of  (76)  are  both  .955. 


Evaluation  of  (12)  gives 

u,=  i,  v,= 


.wm 

.2W4 
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These  matrices  and  their  traces  and  eigenvalues  are  unequal,  but  their 
determinants  are  both  1. 
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3.  UNKNOWN  CORRELATION 


In  this  section,  the  correlation  values  {RjJ  are  unknown,  and  the  only 
information  available  about  the  random  process  is  a finite  set  of  N data 
points  , X2,  ...,  X^,  from  which  we  have  removed  the  sample  mean.  From 
these  N data  po./its,  we  desire  an  estimate  of  the  spectrum  G(f).  But  we 
cannot  minimize  or  utilize  any  ensemble  averages  as  was  done  in  section  2, 
since  we  have  only  a finite  segment  of  one  member  function  to  work  with. 

3.1  PHILOSOPHY  OF  APPROACH 


For  the  known  correlation  case  above,  we  had  the  set  of  normal  equations 

-O  <7w> 

I , p> 

£b!X„  - 1L  (786) 

W*  | J 

where  and  {bH;  were  the  unknowns.  Now  in  the  unknown  correlation 

case,  we  make  a change  by  assuming  that  and  Byf  are  known*  (along  with  "Km 

for  1*1  < p-l , from  lower  order  solutions),  and  by  letting  l?f  andl?^  be  unknown. 
The  equations  in  the  unknowns  are  still  linear,  and  the  solution  is  given  by 


/?-A v-W 


} I — ^ — J>“l 


(pi2)  . 


(79A) 


(79B) 


*The  manner  of 


specifying  and 


will  be  considered  in  subsection  3.4. 
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(80A) 

V ^xf  . 

(80B) 

(It  must  be  noted  that  in  this  section  denotes  an  estimate  of  the  true 
(unknown)  correlation  value;  for  notational  convenience,  no  distinguishing 
symbol  has  been  added  to  to  emphasize  this  distinction.)  However,  we 
shall  insist  that  the  correlation  estimates  (80)  that  we  obtain  at  the  p-th 
stage  satisfy 


(81) 


in  keeping  with  property  (2).  Since  equations  (78)  and  (81)  are  identical 
to  those  encountered  in  the  known  correlation  case,  the  mathematical  defini- 
tions and  interrelationships  employed  there  can  be  applied  here  also.  How- 
ever, some  of  the  properties  and  physical  interpretations  may  be  different, 
since  we  are  now  dealing  with  estimates,  rather  than  true  values. 


(82) 


(8j) 


To  solve  (78),  we  begin  by  defining 

k 

Now  consider  p=l  in  (78);  we  have 

fx  - ft  , Btx  - K ■ 

Now  if  and  are  known,  we  can  compute  unknowns  K,  and  R...  But  by 
constraint  (81),  and  jjj0  must  be  chosen  such  that 

AfK  - . 

Thus  when  we  select  ^ and  constraint  (84)  must  be  kept  in  mino;  that  is. 


(84) 
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and  cannot  be  specified  independently  of  each  other. 


At  stage  p(£2),  if  Aj  is  known  (and  \^hr  are  known  from  earlier 
stages  with  property  "R„* 0*k<  p~l ) , we  could  solve  the  linear  equations 
(78A)  for  {$!?}[  ar|d  ftp*  according  to  (79A)  and  (80A),  where  the  lower  order 
quantities  in  (79)  and  (80)  are  available  from  earlier  stages.  Similarly 
if  is  known,  we  use  (79B)  and  (80B)  to  solve  (78B).  However,  by  (81), 

fa)  u 

we  must  constrain  the  selection  of  Al  and  Bl  . 

T P 

To  see  exactly  what  constraint  (81)  implies  about  the  selection  of 
and  , notice  that,  for  p>2,  (and  defining 

- $ (T-  a;b^’)V + 


= % l^p-*  V* 


(35) 


where  we  have  employed  (80A)  and  (79A),  Now  define  forward  extrapolated 
correlation  estimates  based  on  order  p-1  according  to  (see  (52)  and  (53)) 

>->) 


tT.  S.  aVi£!  fc 


W*  I 


(86) 


where 


i 2 t , Os  nsp-l 


(87) 


Then,  in  particular,  the  one-step  forward  extrapolated  correlation  estimate 
based  on  order  p-1  is 

(88) 


a Af’ V = I’AtV  • 


Also  define  (see  (12)) 


\ 
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4 


u 

I; 


V ■ • l®t\ . 


'p-t  ' . (89) 

This  quantity  has  the  physical  interpretation  as  the  estimate  nf  the  corre- 
lation matrix  of  the  (p-l)th  order  backward  residual  at  zero  time  delay 
(see  (30));  its  properties  are  considered  in  subsection  3.3.  Then  by  means 
of  (88)  and  (89),  (85)  can  be  eApressed  as 

4 = Afvr, . <9o> 

(This  equation  is  similar  to  a combination  of  (14)  and  (59)  for  the  known 
correlation  case.) 

At  the  same  time,  by  (80B)  and  (79B)  (and  defining  ^^-l), 


= SbrV 

*»««  > r ;-o  o -j 


(91) 


Itow  define  backward  extrapolated  correlation  estimates  based  on  order  p-1 
as  (see  (55)  and  (56)) 


where 


Tt'.SBrc  t 

, OS.s  p-l 


(92) 


(93) 


Then,  in  particular, 
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siT-lc^lrv 


(94) 


Also  define  (see  (12)) 


'-v.  ■ • I*'*. 


(95) 


This  quantity  is  an  estimate  of  the  correlation  matrix  of  the  (p-l)th  order 
forward  residual  at  zero  time  delay  (see  (30)).  Then  by  means  of  (94)  and 
(95),  (91)  can  be  expressed  as 


(This  equation  is  similar  to  a combination  of  (14)  and  (60)  for  the  known 
correlation  case.)  But  now  it  can  be  shown  (see  appendix  G)  that  the 
extrapolated  correlation  estimates  in  (8b)  and  (94)  satisfy 

Therefore,  if  (81)  is  to  be  satisfied,  (90)  and  (96)  in  conjunction 
with  (97)  force 


(96) 


(97) 


(98) 


(This  reduces  tc  (84)  for  p=l.)  Thus  the  selection  procedure  of  and  ^ 
at  the  p-th  stage  must  be  done  according  to  (98),  where  and  are 
quantities  already  available  from  the  (p-l)th  stage,  according  to  (89)  and 
(95).  The  precise  selection  procedure  will  be  undertaken  in  subsection  3.4. 
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3.2  COMPARATIVE  FEATURES 


There  are  alternative  techniques  to  the  estimation  of  the  correlation 
matrices  and  the  spectral  density  matrix  that  could  be  considered.  For 
example,  the  standard  Yu’e-Walker  technique  (e.g.,  Ref  2,  page  186)  uses 
correlation  estimates 


aT  1r^  ) 


(99) 


where  the  sum  is  over  all  nonzero  summands,  and  then  solves  recursively  for 
m:  and  ffl  via  the  method  in  subsection  2.1.  This  apriori  decision 
on  the  form  (99)  of  the  correlation  estimate  gives  poorar  spectral  estimates 
for  M=1  (Refs.  2 and  3),  and  probably  does  so  for  M>1.  The  estimated 
correlation  matrix  is  Hermitian,  block  Toeplitz,  and  nonnegative 

definite: 


where^  isMxl.  However  the  stability  of  the  correlation  recursion  (52) 
is  unknown  to  this  author.  The  estimate  (99)  is  unchanged  by  the  addition 
of  more  stages,  that  is,  larger  values  of  p. 
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Another  technique  would  be  to  minimize  the  prediction  error 

^ * ^A„X  p+i  ski  M 

p 

over  the  available  data  points  directly,  by  choice  of  • We  have  the 
error  matrix 

r^< Y*  Y*H  = AwH 

I 1 


(101) 


(102) 


where 


_J — Y y H 

K^l  Ah‘’’  A 


0 * w»;  r?  < f. 


(103) 


The  optimum  coefficients  for  minimum  trace  of  the  error  matrix,  (102),  are 
solutions  of 


±/fc- 


Otrt 


(104) 


Matrix  [5mnJ?  is  not  block  Toeplitz,  and  a significant  computer  problem 
exists  for  M>1  when  it  is  noted  that  solution  of  linear  equations  (104)  must 
be  done  anew  for  each  different  value  of  p.  This  was  a good  technique  for 
spectral  estimation  when  M=1  (see  Ref.  3);  however,  computer  time  was 
greater  than  for  the  Burg  technique.  Moreover,  stability  of  the  correlation 
recursion  (52)  is  unlikely  in  view  of  the  (occasionally  unstable)  results 
for  M=1  in  Ref.  3. 

This  technique  could  be  extended  to  include  backward  prediction  in 
addition  to  (101).  However,  the  lack  of  the  block  Toeplitz  property  and 
lack  of  stability  make  it  a very  undesirable  technique. 


\ 
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The  technique  suggested  here  (in  subsection  3.1)  lets  the  correlation 
estimate  be  yielded  according  to  solution  (80),  once  partial  correlation 
coefficients  and  have  been  specified.  And  we  shall  see  in  subsection 
3.4  that  these  latter  quantities  are  determined  according  to  a physically 
meaningful  minimization  problem.  Stability  of  the  correlation  recursion 
(52)  has  not  been  proved;  however,  numerous  examples  have  all  yielded  stable 
solutions.  The  estimate  (80)  is  unchanged  by  the  addition  of  more  stages, 
that  is,  larger  values  of  p.  And  it  will  be  seen  that  the  current  technique 
reduces  to  Burg's  algorithm  (Ref.  4)  for  M=1 . Thus  the  current  technique 
appears  to  be  very  attractive  among  those  techniques  that  employ  an  all -pole 
representation  of  the  input  process. 


3.3  PROPERTIES  AND  INTERPRETATIONS 


The  quantities  U^.(  and  \£_(  were  defined  in  (95)  and  (89)  and  were 
interpreted  as  estimates  of  the  correlation  matrices  of  the  (p-l)th  order 
forward  and  backward  residuals,  respectively,  at  zero  time  delay.  It  is 
shown  in  appendix  H that  they  satisfy  the  recurrence  relations 


, u.-*, 

Vf*(r-^ApVr , V.--K  , 


(105) 


just  as  for  the  known  correlation  case.  It  is  also  shown  that 


and 


HH>  %H-vr> 


Up  - W Vp  . 


(106) 


(107) 
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However,  we  are  not  able  to  prove  or  Vp  nonnegative  definite  without 
specifying  the  method  by  which  and  are  selected;  no  relations  like 
(30)  and  (32)  exist  here. 


By  means  of  (106),  the  constraint  (98)  on  selection  of  and  takes 


the  form  (see  42)) 


A*Vur 


(10£ 


This  will  be  used  in  the  next  subsection. 


3.4  EVALUATION  OF  PARTIAL  CORRELATION  COEFFICIENTS 

We  recall  from  subsection  2.1  that,  in  the  known  correlation  case,  the 
partial  correlation  coefficients  and  gjf  minimized 


tr  «»d 


(109 


respectively,  when  lower  order  stages  had  already  been  optimized.  We  extend 
this  idea  to  the  unknown  correlation  case  as  follows:  let  (as  in  (18)  and 


(19)) 

'C  * x„ , X,  , i-  M nL 


and  for  p*  I , define  errors  (resiouals) 


I r 


i * k £ N 


(11 


The  block  diagram  for  (111)  is  identical  to  that  in  figure  1 on  page  7. 
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Define  for  psl,  the  error  (residual)  matrix  over  the  available  data 
points  as 


fcP  N-J>  'K  <K  *7  ' 


this  nonnegative  definite  matrix  is  a.i  unbiased  estimator  of  Y* 
Substitution  of  (111)  in  (112)  yields 


where 


Substitution  of  (111)  in  (115)  yields 

. (P)  Ini)  IrtH 


(112) 


_ cT  -v  A'fiH 
rip  Op- 1 «Op-»  Hp  r 

Af’cr^  A^H 
'Y  3r,  , 

(113) 

4™  - 

•-'p-i 

_j ^ _ 

N-f  K^h  ’*  'h 

«v  War 

°p-.  , 

(114A) 

oS^  » 

op-1  “ 

_i . Jr  yr‘)7<r^ 

N-p  CsfTl  ‘k  A".  ) 

(114B) 

in 

i Jr  *3H  7V*‘>h  _ 
N-f  jT^i  Am  Ac-. 

"V*  * 

(114C) 

pal  , error  matrix 

= — i — J:  2'p)2ph  . 

(115) 

Now  error  matrices  Ep  and  are  Hermitian  and  nonnegative  definite. 
Therefore  matrix  M-A.  is  Hermitian  and  nonnegative  definite  for  any 
MxM  weighting  matrix  A0: 


(116) 


(117) 


•>! 

Vf 

s 

-'■J 
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« 

1 

? 


for  any  Mxl  matrix  V*  Also  since 

4r(AaEpAH)  r (118) 

only  the  product  AA^Ax  matters  in  so  far  as  the  trace  of  \EfJ^  is  concerned; 
notice  thatA  is  Horan  ti an  and  nonnegative  definite.  We  shall  be  interested 

H K 

in  minimizing  the  traces  of  weighted  error  matrices  A,Epjta  and  1 f>  ij  ; 
the  exact  choice  of,  and  the  reason  for,  weightings  and  Q will  be  under- 
taken in  the  next  subsection. 


i \ tp) 

Now  if  we  were  to  minimize  £r(Ap.,E^  by  choice  of  Ap  , we  would  find 
(see  appendix  B for  method)  that  we  must  solve 


(119) 


and  the  choice  of  Ap_,  would  be  irrelevant.  Also,  if  we  were  to  minimize 
Mi:sF)  by  choice  cf  Bj?,  we  would  find  that  we  must  solve 


r 

1 p-,  £)p 


n 

' P-i 


(120) 


and  the  choice  of  would  be  irrelevant.  Furthermore,  we  would  not  satisfy 
constraint  (108)  generally.  But  since  the  behavior  of  error  matrix  Ff  is 
just  as  important  as  that  of  Ep,  we  should  take  both  matrices  into  account 
in  any  erre"  minimization;  in  fact,  for  known  correlation,  recall  that  the 
determinants  of  residual  matrices  and  were  equal. 


I 


J , 


\0 

\ 

jr< 


r < 


6’ 


e i 

t- 


.o 

■'i 
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We  therefore  choose  to  minimize  the  sum  of  the  traces  of  the  weighted 
error  matrices 

(,2,) 


where  and  are  Hermitian  and  nonnegative  definite, 

by  choice  of  and  subject  to  constraint  (108).  If  we  let 

f'X  - % < «• , 

then  we  can  express 


(122) 


Ej  +■  Fr,  ^ = 

‘r  Vf.,  [sjj?  - G*  Df.'  - 5*f up  ftf  4 ffr  ur,  ej 

in  terms  of  the  single  unknown  matrix  Gp.  Our  problem  therefore  is  to 
minimize  the  trace  of  (123)  by  choice  of  the  single  quantity  Gp,  subject 
to  no  constraints;  we  can  then  solve  for  the  best  coefficients  according  to 


f-spv;; , 


(124) 


Also  we  can  compute  the  correlation  estimate  from  (90)  and  (88)  according 

\ - *r)+  ■ 


(125) 


In  appendix  I,  it  is  shown  that  the  minimum  of  the  trace  of  (123)  is 
realized  when  Gp  is  the  solution  of  the  bilinear  matrix  equation  (Ref.  9) 

ffj,  f ^ Gf  ~ ^ +? ; 


(126) 
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where 


*-V'#£rr; 

f - V'  ^ ur:  < 

* - a;:  v£  s^f  . 

Uniqueness  of  the  solution  uf  (126)  is  considered  in  subsection  3.6.  (It 
is  interesting  to  note  that  the  separate  minimizations  in  (119)  and  (120) 


yield 


Gp°<-  y - v~  - o • 


Thus  whereas  both  these  quantities  had  to  be  equal  separately  to  the  zero 
matrix,  we  now  require  only  that  they  be  equal  to  each  other.) 


For  the  special  case  of  M=!  (a  univariate  process),  (105)  and  (108) 


yield 


¥ _ aV1* 


(m=  I). 


Then  (126)  and  (127)  can  be  solved  for  the  scalar 


ft . u 

I JM  ' Jlp-I 


(m=  ') 
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Now,  if  and  only  if 

(m  = j),  (i3i) 

(130)  reduces  to  Burg's  algorithm  (Ref,  4);  in  fact,  it  can  be  shown  that 

(131)  is  the  only  choice  of  weights  in  (130)  which  guarantees  a stable 
correlation  recursion  for  M=l.  Thus  we  shall  insist  that  the  weights 
satisfy  (131)  when  we  deal  with  their  selection  below. 


3.5  WEIGHTING  OF  ERROR  MATRICES 


It  is  necessary  to  apply  weighting  to  error  matrices  Ep  and  Fp  in  (112) 
and  (115),  prior  to  minimization  of  the  trace  in  (121),  for  several  reasons. 
First,  without  weighting,  the  larger  amplitude  components  of  errors  (111) 
would  receive  most  of  the  emphasis  in  the  minimization;  thus,  some  weighting 
inversely  proportional  to  the  component  strengths  is  desired.  Second,  it 
is  desired  that  stable  correlation  recursions  result  and  that  matrices  Up 
and  Vp  be  nonnegative  definite.  Without  weighting,  it  has  been  discovered 
(by  an  example  to  be  presented  in  subsection  3.9)  that  both  of  these  require- 
ments can  be  violated.  Third,  we  will  insist  that  the  scaling  property 
introduced  in  subsection  2.2  hold  for  the  unknown  correlation  case  as  well; 
that  is,  if 


(132) 


we  shall  insist  that  the  coefficients  satisfy 


a:=j>a!V] 

Iff  » D Bt’D"' 


(133) 
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The  matrix  equation  (126)  can  be  combined  with  (122)  to  yield  the 
simultaneous  set  of  equations 


(134) 


We  now  consider  several  possible  choices  of  weightings  and  f^_.  that  tend 
to  simplify  the  form  of  (134).  The  first  choice  is  no  weighting: 


— > fp-j  _ X . Cl»o»ce  I (135) 

The  problem  with  this  choice  is  that  the  weighting  is  not  related  to  the 
error  component  strengths,  and  it  may  be  readily  verified  that  the  solutions 
to  (134)  and  (135)  do  not  satisfy  the  scaling  property  (133).  Also  an  unstable 
correlation  recursion  can  occur.  However,  the  solutions  do  reduce  to  Burg's 
algorithm  for  M=l;  see  (131). 


Our  next  candidate  weighting  is 

n \ i "I 


V ■ H-' , rP-,  - V 


Choice  2 


(136) 


which  are  Hermitian  and  are  nonnegative  definite  if  Up_-j  and  Vp  ^ are 
nonnegative  definite.  This  weighting  is  inversely  proportional  to  the 
component  strengths,  as  desired;  more  will  oe  said  on  this  below.  The 
equations  (134)  become 


* s£e f . 2 s, 

o- 


f>“‘  > 


(137) 


The  solutions  of  (137)  satisfy  the  scaling  property  (133),  and  they  reduce 
to  Burg's  algorithm  for  M=l;  (129)  shows  that  (131)  is  satisfied  for  the 
choice  (136).  Although  stability  of  the  correlation  recursions  (52) 
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and  (55),  and  nonnegative  definiteness  of  Up  and  Vp,  have  not  been  proven 
for  genera]  M£2,  no  counter  examples  have  been  discovered. 


We  next  consider 


r * v 

T< 


% 


-i 


Clio>Ve  d 038) 


in  which  case  (134)  becomes 

A?  + B f-  5 

o. 


(139) 


However,  the  weighting  (138)  is  not  necessarily  Hermitian,  is  not  necessarily 
nonnegative  definite,  and  is  not  directly  related  to  the  error  component 
strengths.  Also  the  solutions  of  (139)  do  not  satisfy  the  scaling  property. 
Furthermore,  the  solutions  do  not  reduce  to  Burg's  algorithm  for  M=l,  and 
can  yield  unstable  correlation  recursions  for  M=l. 


The  last  choice  is 


Ctal 


ce 


(140) 


which  are  Hermitian  and  nonneqative  definite,  and  for  which  (134)  becomes 


Ap  Vp.,  i-  By  - 5^5^  Vp.,+  jp.,  ^ , n41) 

A^’V  Up.,  B)f  - o. 

This  choice  is  a very  interesting  one  in  that  the  solutions  of  (141)  are 
immediate  and  do  not  require  that  a bilinear  iiiu  tr  i x equation  be  solved.  The 
weighting  (140)  is  inversely  proportional  to  the  error  component  strengths, 
and  the  solutions  of  (141)  d£  satisfy  the  scaling  property.  In  fact,  this 
choice  is  very  close  to  Choice  2,  since  Up and  S^are  both  estimates  of 
the  correlation  matrix  of  process  at  zero  time  delay,  and  should  be 
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fairly  close  to  each  other.  However,  the  solutions  of  (141)  do  not  reduce 
to  Burg's  algorithm  for  M=l,  and  the  correlation  recursion  (52)  can  be 
unstable,  even  for  M=1 , In  fact,  the  solutions  to  (141)  are  identical 
to  those  for  Choice  3 for  M=1 . 


Therefore,  of  the  four  choices  considered,  only  Choice  2 in  (136) 
yields  solutions  that  satisfy  the  scaling  property  (133)  and  reduces  to 
Burg's  algorithm  for  M=1 . The  stability  of  the  correlation  recursions  has 
not  been  proved  or  disproved  for  choice  (136)  of  weighting. 

There  is  another  strong  reason  for  choosing  weighting  (136),  which  has 
to  do  with  a whitening  interpretation.  We  recall  that  U^-(  and  \/  } 

defined  in  (95)  and  (89),  are  estimates  of  the  correlation  matrices  of 
processes  fY**0}  and  » respectively,  at  zero  time  delay.  Now  let 

(for  non-negative  definite  and  V?-,  ) 

u vl”  v =v  vh  ;u: 

vp-’  v7-»  j V*  V'  V*  > 

where  and  V^,  are  (lower  triangular)  square  root  matrices.  Ther.  scaled 
processes 

, ir-v;r.  „« 

each  have  estimated  correlation  matrices  at  zero  time  delay  equal  to  I;  that 
is,  all  the  components  of  (or  ) have  unit  power  and  are 

uncorrelated  with  each  other  at  zero  time  delay. 


Now  def i ne . for  p+ 1 i )c  < N , 


f ■ «;:  • <j;  hr-  <0  ■ tr-  k sr , 
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where 


(145) 


Also  define  the  estimated  correlation  matrix  at  zero  time  delay  of  process 

^ ' W ' V ■ ■4-’’ H-'  > 


(146) 


where  we  have  used  (144)  and  (112).  Therefore 

H - ^(v;:\-;£.)=  ^e,),  (wi 

— I 

where  we  have  used  (1-1)  and  (142).  Thus,  minimizing  the  trace  of  iy,  e,. 
by  choice  of  A^  ,is  equivalent  to  minimizing  the  trace  of  fp  by  choice  of 
^ (see  (144)),  where  process  fy^is  the  error  in  prediction  of  (p-l)th 
order  processes  with  estimated  correlation  matrices  at  zero  time  delay  equal 
to  I. 


In  a similar  fashion,  for  pt-Kfcs  N, 


zT-ifvr) 


(148) 


where 


Ef ' . V,::  wr-.  . 


(149) 


And 


% ■ if&W- v;f,  vr: 


(150) 


with 


hr  ^ - "tr  ( Vp_,  Fy). 


(151) 


36 


TR  5501 


If  we  solve 
(108)  along  with 


(145)  and  (149)  for  and  Bp*  , and  then  utilize  constraint 
(142),  we  find  that  the  constraint  takes  the  form 


(152) 


This  could  be  used  as  the  starting  point  in  a minimization  of  error  matrices 
and  $ . In  fact,  if  we  minimize  the  unweighted  trace  offp*3p  by  choice 
of  hf  , we  *ind  the  optimum  choice  to  be  given  by 


where  the  notation  is  an  obvious  modification  of  (114).  By  employing  (145), 
(143),  and  (142),  we  can  show  that  (153)  is  equivalent  to  (137),  as  it  must 
be.  (This  alternative  approach  may  be  useful  for  proving  the  stability  of 
the  correlation  recursion.) 

3.6  SOLUTION  OF  BILINEAR  MATRIX  EQUATION 

If  we  substitute  definitions  (127)  into  bilinear  matrix  equation  (126), 

*/»  */j 

and  premultiply  by  JL(  and  postmultiply  by  , we  obtain  the  equation 


Gj,  * + 2 §r  = f*  * <,54> 

where 

r : p*1 

- “V*  7 7- 

* ' W 4?  v;  r* 

(155) 

p . -V  if'  sg  rft 
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* i 


Now  the  Hermitian  matrices  v and  jj  are  non-negative  definite;  e.g., 

1 M - (v; rrTv)h^(\:  [;» > o <« 

for  any  MX1  matrix  ^/,  since  £ is  non-negative  definite.  We  have  employed 
the  Hermitian  property  of  \^_t  and  j^,  above;  see  (118)  et  seq.  This  means 
that  the  eigenvalues  of  % and  j?  must  be  non-negative.  Therefore  the 
solution  of  (154)  exists  and  is  unique  (Ref.  10,  eq.  3). 


(156) 


Solution  of  the  bilinear  matrix  equation  (126)  or  (154)  has  been  addressed 
by  many  authors  (Refs  9 - 17).  In  particular,  for  the  equation  involving  MxM 


matrices, 


X B+  A X = C, 


(167) 


one  form  for  the  solution  is  given  by 


X = PQ  , 

cjt, 

are  MxM  matrices.  The  constants  are  given  by  (Ref.  18,  pp.  87-88) 


and  the  matrices  are  given  by 


Ah-AA^+^x 


(158) 


(159) 


(*osl),  (16°) 


(A  - r)  ■ 


Here,  M-2  full  matrix  multiplications  are  necessary  when  we  note  that  AM  = 0 
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by  the  Cayley-Hamilton  theorem. 


For  M = 2,  (159)  takes  the  form 

T*  c®  - (A-V<o)l)  C 

where  we  have  used  the  Cayley-Hamilton  theorem  to  express 


fer  M=2  , 


(162) 


8*  ^ -trig)  £ - detfe)  X for  h--2 


(163) 


Equations  (162)  and  (158)  are  the  forms  used  in  the  FORTRAN  program  for  M = 2. 
3.7  SPECTRAL  ESTIMATION 

Having  obtained  correlation  estimates  by  means  of  (82)  and  (80A), 

we  now  extrapolate  these,  as  in  subsection  2.3  (equations  (52)  and  (66)),  to 

yield  V . y 

rv*  K^_„  ? pt)  i 

(164) 


w-  l 


^ = R-n,  } m v 0 . 

This  defines  an  infinite  sequence  which  is  assumed  stable;  its  Fourier 

transform  will  be  taken  as  the  spectral  estimate  of  the  process  under  consid- 
eration. In  a manner  identical  to  that  given  in  appendix  F,  it  is  found  that 


J_  (165) 
2A  ' 


= aiffff'u,  if) 

n)2  *0  j ^ ) 

where  Op  and  H^ff)  are  given  by  (95)  and  (68),  respectively.  It  follows  that 

j. 

2* 

r # « w ■ k - Sample  poiotr  (So).  (166) 


Also,  as  in  subsection  2.4,  an  alternative  factorization  is  available  as 

ifW  - -H 


V i£’(f)’  , If) 


(167) 
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M . 

where  Vj>  and  Hp(f)  are  given  in  (89)  and  (70).  If  or  is  non-negative 
definite,  then  Cr^ff)  is  non-negative  definite,  as  desired  for  a spectral 
estimate.  Since  (165)  and  (167)  are  equal,  we  concentrate  henceforth  on 
form  (165). 

Since 


(168) 


(165)  can  be  expressed  <ij 


(169) 


Since  fi^jFlis  Herrnitian,  matrix  Cr^tf)  need  be  computed  only  on  ana  above  its 
main  diagonal,  at  each  frequency  of  interest.  Efficient  computation  of 
Hjfff)  by  means  of  an  FFT  is  undertaken  in  appendix  J.  It  is  shown  that  we 
need  to  perform  M*  fy-point  FFTs  of  p+1  non-zero  numbers,  in  order  to  evaluate 
H^(f)  at  Np  frequency  cells  in  the  frequency  range  ^ 


Real  Multivariate  Process 

The  results  above  have  been  derived  for  a complex  multivariate  process 

..  jp 

[XK]  • For  a real  multivariate  process,  Ij,  is  real  and  [A^}d  are  real  • Then 

$(-f)  = whf  for  a real  process,  (170) 


and 


for  a real  process. 


(171) 
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Thus  we  need  compute  matrix  6r^(f)  only  for  f > c> , for  a real  multivariate 
process. 


In  order  to  avoid  complex  matrix  multiplications,  we  develop  (169)  more 
explicitly;  let 


H = \ If )+.  !,«, 


(172) 


where  J^(f)  and  l^(f)  are  real  MxM  matrices  at  each  f.  Then  since  U^,  is  real 
\Jp  , and  upon  suostituting  (172)  in  (169),  we  find 


where 


fcfHvAwW*  .•  Mif)- 1 «ify] 

Q rto|  yrtctn ; 


(173) 


M(f)-*l(f)l^K^)T  (174, 

Since  M(f)  is  real,  the  quantity  (M(f)-i  M(f)  is  zero  on  the  main  diagonal; 
therefore  we  need  not  compute  the  main  diagonal  of  M(f).  All  the  matrix 
multiplications  in  (173)  are  real. 


Real  Bivariate  Process 

We  now  further  specialize  to  M = 2,  a bivariate  process.  Let  the  real 
and  imaginary  parts  of  the  filter  transfer  function  be  denoted  by  XX 
and  YY,  respectively  (where  these  symbols  are  unrelated  to  X and  Y introduced 
earlier);  that  is 

HT(F)»  XX(f)+.-YV(t).  "75> 

Then  from  (172),  for  2X2  matrices, 

\(i)  ^tfAdjUKV)]-  U fcfofBj*  H XX(f) . X\(f),  (He, 
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l^fr) ' Lt  \ftfj  ^ (f)J  -Adj  I^^ff^j  = Aclj  Wtf) E (f) . 

Substitution  of  (176)  and  (177)  in  (173)  yields  spectral  estimate 

<?^  • H$f)f  [xx^ij,  '^Wt+^  (f) n(f)T+ » Hf)-  i WT] 

■&r  a r»a!  ^ivanWe  process. 


077) 


where 


Mlf)*YXW  ‘4  »dff)r 


(178) 


(179) 


The  2X2  matrices  involved  in  (178/  are  all  real,  and  XX/\(f)  and  YY.(f)  are 
the  adjoints  of  the  real  and  imaginary  parts  of  , respectively.  The 

form  (178)  is  used  in  the  program  for  the  spectral  estimate  of  a real 
bivariate  process. 

3.8  TERMINATION  PROCEDURE 

For  unknown  correlation,  the  correct  value  of  p to  use  in  (79)  and  (80) 
is  unknown.  We  adopt  the  Akaike  information  criterion  (AIC)  derived  in 


Ref.  19,  page  719: 


A 4.  Cp  = N ■ h 4*4  \jp  f 2 *1  j> 

aN  *n  <ki  \ f 2 h7f>. 


(180) 


where  we  have  utilized  (107);  namely,  we  compute  AICj  for  p = 0,  l,..,p  , 

K max 

ana  we  use  that  value  of  p,  pL  ^ , for  which  AIC  is  a minimum.  Selection 

best  p 

of  p is  discussed  below, 
max 

For  purposes  of  updating  Up  and  V^,  we  can  combine  (105),  (106),  and 
(122)  to  yield 


- . vAc 

\V 
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Up  = Uf_,  - A,  Sj,M , Vp  = vp_,  - Bp  Gj, , ('si) 

in  terms  of  the  solution,  Gp,  of  bilinear  matrix  equation  (126). 

At  this  point,  it  is  worthwhile  to  review  the  procedure  adopted  here. 

From  the  actual  data,  we  could  have  estimated  the  input  correlation  matrix 
via  (99)  (or  some  scaled  version  of  it).  Also  we  could  have  used  (112)  and 
(115)  as  error  matrix  estimates;  in  fact,  these  matrices  are  guaranteed 
Hermitian  and  non-negative  definite.  However,  since  <Jet  Ej,  1 det  we  would 
have  had  to  settle  on  some  average  like 

I"  ( * W lj~)  - £ ( In  cW  Ej,  + k det  (182) 

for  purposes  of  che  information  criterion.  As  for  the  spectral  estimate,  we 

I H 

could  have  adopted,  instead  of  (165),  the  quantity  A^ff)  Ep  ' , 

Fp  Hg  (f)  , for  example. 

Instead,  we  have  chosen  consistently  to  stick  with  the  results  of  the 
normal  equations  (78).  Thus  the  estimate  of  the  input  correlation  matrix  is 
obtained  from  (80)(and  (82));  the  estimates  of  the  correlation  matrices  of  the 
residuals  are  given  by  (89)  and  (95)  (or  more  computationally  convenient  yia 
(181));  and  the  spectral  estimate  is  given  in  terms  of  or  by  (165)  or 
(167),  respectively,  for  P=Pbest-  The  major  gap  in  this  procedure  is  that 
we  have  not  proved  that  Up  or  is  non-negative  definite  for  Choice  2 of 
weighting  in  (136);  however,  no  counter  examples  have  been  discovered. 

Our  selection  of  Pmax  ’s  accomplished  as  follows:  in  ref.  1,  page  575, 

Akaike  is  quoted  as  suggesting  Tor  M = 1,  a univariate  process.  Since 

the  number  of  coefficients  evaluated  is  p,  and  the  number  of  available  data 
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points  is  N,  this  ratio  was  upper  bounded  by  ^le  extend  this  idea 

directly  to  the  multivariate  case:  the  number  of  scalar  coefficients 

evaluated  is  , and  the  number  of  available  scalar  data  points  is  MN. 
Upper  bounding  this  ratio  by  3N  , we  find  we  should  choose  the  filter 
order 


M (183) 

in  terms  of  the  number  of  data  points,  N,  and  the  dimensionality  of  the 
time  series,  M. 


3.9  EXAMPLES 

It  is  worthwhile  to  summarize  here  the  sequence  of  calculations  required. 

For  data  X-| , X^,..,  available  (with  the  sample  mean  removed),  we  have 

\ 

'K 


nC’2T-  X.,  '-'**''*  H 

Mo  i <"  V 
5*  = N-l  kH  **  **  " 0 


Then  forp>i  and  choice  (136)  of  weighting, 

V,  - > w 
■<  * Vp-i  Sp-, 

f - u,:: 

r -- » * sfi 

(r.  via  (l  ifc) 


\i,  ■ Vp.,- % 

ATCpvAi 


(184) 


085) 


h 

5! 

n 

ft 

u 


1 


L 

V 


1*1 
* 4 

i. 

't 


?, 

j 

\§ 

; o 

* < 
* v 


l 

y 


& 

c 
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j 


I -*jr  v4pl  v^n 
°r  = w-r i J$5  '*  >*  * 


I jr  3*»  7,»)H  a J___  r:  7»y>M=  <$»' 

N-p-i  **''  ^-'  /J-f-l  Jcr*2*  '~7 

I ^ V^W* 

N_f- i fc.^j  'k  £*-<  • 


086) 


(187) 


For  p = pmax’  it:  is  not  necessary  t0  compute  (186)  through  (187).  uu.c.i  the 
je  oi 

estimate  (165). 


best  value  of  p,  p is  found  from  AIC  , we  can  then  compute  the  spectral 
best  p K 


We  now  consider  an  example  for  M = 2,  N = 4: 


A* 


-i 

-2 


A- 


4,  - V,  = K„  = 


Then  for  weighting  (136),  we  find 


I 0 IS 
I.S  2 s 


(188) 


n i? 
if  if 


aV>  = _ j_ 
; ft,  = O.  = 6 


I J 
( 1 


(189) 


The  eigenvalues  of  are  (-3tfT})/i2  , which  are  both  bounded  by  1 in  magnitude, 
as  they  must  be  for  the  correlation  recursion  (164)  for  n = ] to  be  stable. 
Also, 

r~  - 

6 


G 


IQ 


(190) 
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which  is  non-negative  definite.  Thus,  for  weighting  (136),  all  the 
desirable  properties  are  realized. 


However,  for  no  weighting,  (135),  we  find,  for  the  same  example  (188), 


\to  30]  m w ± 

J>  V®.’  + 


-10 

-V 


0 


(191) 


The  eigenvalues  of  A?  are  4/9  and  -10/9;  since  the  latter  is  larger  than  1 


in  magnitude,  the  recursion  , n?l,  is  unstable.  Also 

v, 


[” 


(192) 


which  is  not  a non-negative  definite  matrix.  It  is  found  that  the  spectral 
estimate  obtained  from  (165)  has  frequency  rangejwhere  the  two  autospectra 
(diagonal  terms  of  (165))  are  negative,  and  where  the  magnitude-squared 
coherence  can  be  negative  or  greater  than  1.  These  are  all  unacceptable. 


For  the  alternative  example  for  M = 2,  N = 4,  of 

vl;; 

and  no  weighting,  we  find  a stable  correlation  recursion,  but  and  are 
not  non-negative  definite,  and  values  of  the  magnitude-squared  coherence 
greater  than  1 are  realized  in  some  frequency  ranges.  Because  of  these 
unacceptable  behaviors,  the  choice  of  no  weighting,  (135).,  is  discarded  from 
future  consideration. 


, x,* 


- 

-l.is 

si 


.75 
M 3 


(193) 


An  example  for  M = 2,  N = 100,  and  weighting  (136),  generated  via  (71)  - 

(73)  of  subsection  2.5  yielded  the  results  below;  the  program  and  its  output 

are  given  in  appendix  K.  We  find  p =1  and 

best 
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<= 


.H/5! 

-11024 

' f. 

• 5£4>I3 

.770'** 

.63432 

o603s| 

■ ' I 

- 63142 

. U 573 

(194) 


U,  = .09110 


°n(,U  - C0f47 

■corn  i.  02  364 


1*111  .3W89~1 

V =.  OlHo  I 

' [mn  ).  ioo#7j  ♦ 


(195) 


It  is  worthwhile  to  compare  these  estimates  for  N = 100  with  the  exact  values 
in  (76)  and  (77).  The  scale  factor  .09110  in  (195)  is  unimportant  and  is  due 
to  the  fact  that  the  white  noise  used  here  had  variance  1/12  rather  than  1 
as  in  (73);  except  for  the  scale  factor,  the  matrices  in  (195)  have 
determinants  equal  to  1.  The  estimated  magnitude-squared  coherence  reaches 
a maximum  of  .999745,  versus  the  true  peak  of  .999013. 

Observations  from  other  examples  of  real  bivariate  processes  have  pointed 
out  that:  the  eigenvalues  of  A,  and  8*  are  identical  and  are  bounded 

by  1 in  magnitude;  the  eigenvalues  of  and  are  not  identical  for 

. (rt 

p>2  . and  can  be  larger  than  1 in  magnitude;  and  the  eigenvalues  of 
and  bT  for  n<p  can  be  larger  than  1 in  magnitude. 

Timing  Results 

Some  sample  execution  times  on  a UNI  VAC  1108  for  SUBROUTINE  PCC, 
which  evaluates  the  partial  correlation  coe*  icients,  are  presented  below 
for  M = 2,  a bivariate  real  process. 
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Table  1.  Timing  of  Subroutine  PCC 


N 

pmax 

Time  of  Execution  (sec) 

100 

10 

0.25 

100 

15 

0.35 

1000 

10 

2.63 

1000 

40 

9.23 

10000 

50 

120 

10000 

150 

326 

The  execution  time  is  almost  linearly  propori/onal  to  N and  Pmax.  The 
execution  time  for  PEFTF  was  1.25  seconds,  and  that  for  SDM  was  0.55 
seconds,  both  for  Np  = 1024  frequency  cells;  see  appendix  K for  program. 

4.  SUMMARY 

A method  for  multivariate  linear  predictive  spectral  analysis, 
employing  weighted  forward  and  backward  averaging,  has  been  presented  and 
programmed  in  FORTRAN.  The  method  constitutes  a generalization  of  Burg's 
univariate  algorithm  (Ref.  4)  to  the  multivariate  case. 

The  choice  of  weighting  in  the  error  minimization  is  very  important, 
and  several  candidates  have  been  considered.  The  weighting  retained,  (136), 
is  the  only  one  of  those  considered  that  satisfies  both  the  scaling  property 
(133)  for  all  M,  and  reduces  to  Burg's  algorithm  for  M = 1.  Also,  the 
weighting  retained  is  equivalent  to  minimizing  the  unweighted  traces  of 
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error  processes  that  are  the  differences  of  approximately  white  processes; 
in  fact,  (136)  could  be  used  as  the  starting  point  of  the  error  minimization. 

The  major  gaps  in  the  analysis  are  that  we  have  not  proved  that  Up  and 
Vp  are  non-negative  definite,  and  we  have  not  proved  that  correlation 
recursion  (164)  is  stable;  however,  no  counterexamples  have  been  encountered. 
The  major  analytical  block  in  this  endeavor  is  the  bilinear  matrix  equation, 
(126),  which  requires  special  treatment  for  its  solution. 
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Appendix  A 


PROPERTIES  OF  A SPECTRAL  DENSITY  MATRIX 


Suppose  an  arbitrary  linear  filter  with  impulse  response  {Hn>  is 
excited  by  input  {X^}.  The  output  at  time  kA  is 

where  the  sum  is  over  all  non-zero  summands.  and  Y^  are  M x 1 matrices, 
whereas  Hn  is  M x M.  In  steady  state,  the  spectra  of  the  processes  in  (A-l) 


(A-l) 


are  related  by 


frTW  = Htf) 


(A-2) 


where  transfer  function 


Htf)  * ir  e*p(-» , 


(A-3) 


and  f frequency  in  Hz  and  is  real. 
Now 


£ (f)H  • a • AXe^p(-i2n-f»»A)Rn  = 


(A-4) 


where  we  have  employed  (2).  Thus  ( f ) is  Hermitian  at  any  value  of  f. 
Similarly  Gy(f)  is  Hermitian  at  any  f. 


Also 


C --  X X' 


( A-5 ) 


is  non-negative  definite  for  any  H(f),  because 


- V X V lH^HX  f - 0 w‘6) 


\ 

. 
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for  any  M x 1 column  matrix^  . Therefore 

^ 1 Jf  SY(-f)  * H(f)M  (A-7) 

is  non-negative  definite  for  any  H(f),  l*  then  follows  that 

G*  (f ) is  won  - m eja^ive  de  i i -far  a 1 1 4 . ( A-8 ) 

To  prove  this,  assume  that  Gx ( f i ) is  not  non-negative  definite;  than  if  we 
choose  H(f)*»  I 5 (f -f -j ) , that  is,  an  impulsive  transfer  function  near 
frequency  f ^ , we  get  R^^Gx(f^)  from  (A-7),  which  contradicts  the  conclusion 

M 

that  R must  be  non-negative  definite, 
o 

Thus  a spectral  density  matrix  must  always  be  Hermitian  and  non-negative 
definite  for  all  f.  In  particular,  this  implies  that  all  the  auto  spectra 
(diagonal  terms  of  the  matrix)  must  be  real  and  non-negative.  It  also  implies 
that  all  coherences  are  bounded  by  unity  in  magnitude. 


A-2 
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Appendix  B 


MINIMIZATION  OF  TRACE  OF  ERROR  MATRIX 


From  (4)  and  (5),  we  have 


where 


Let 


X = X-^A„x, 


Q 


\X-c,  XX-Q- 


(B-l ) 


(B-2) 


(B-3) 


Here,£?  is  M x Mp,  \ is  Mpxl , <t  is  MxMp,  and  Q is  MpxMp,  We  notice 

that  Qh  s q.  ana  yMQ°\)-  |1/MXV|1^0  for  any  Mpxl  matrix  VfO,  if  no  exact 

linear  relation  exists  between  the  elements  of  Y , Y : that  is, 

''k-i  ) ’ > 

Q is  Hermitian  and  positive  definite. 


xxH  = (\-a%yx;-'xzaH) 

=i ?.-aeH-eaM+a«Q" 

*K,-£d'c“'(Q-C6")  Qia-caf- 

Let 

U/H] 

fl-CQ"-  / 

%H] 

where Y is  an  Mpxl  matrix.  Then  for  the  M x M matrix  in  (B-5), 

J 


(B-4) 

(B-5) 


(B-6) 


B-l 


tr  <&■  < - j ^ H 

< * * * & ' • * ; • <>- 
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(B-7) 


(B-8) 


The  real  quantity  for  any^?o  , since  Q is  Hermitian  and  positive 

definite;  the  minimum  value  of  Ajj  is  zero  and  is  attained  if  and  only  if 


1/j-O 

Thus  irr\  is  minimized  by  the  choice  of  d as 

J 

-I 


Therefore,  tr  L is  minimized,  attaining  value  zero,  by  the  choice 
H „T7»r 


ea, 


(B-9) 


since  the  leading  two  terms  in  (B-5)  are  independent  of  CL. 


Then  we  have  opt  L = 0 and 

opt  X?  - K-  eQ'e^K-O^c”--  Ko-^0^h  . (B-io) 

Also 

.bfa-ea'cy^-trica'c*)  (B'n) 

It  should  be  noted  that  the  solution  (B-9)  is  attainable  directly  from 
(B-4)  if  the  coefficient  of  & (or(JL)  is  set  equal  to  zero;  this  observation 
will  be  useful  later. 


Equations  (B-9)  and  (B-10)  can  be  developed  as  follows: 


(LlQ,zC  yields»  w1th  ^e  use  (8-2)  and  (B-3), 

R.\ 


R,.. 


K.J 


= &•  Krl> 


(B-12) 


B-2 
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that  is. 


And  (3-10)  can  be  expressed  as 

op* 


’■nr1 


K -ZtiX 

»«  I 


Equations  (B-13)  and  (B-14)  are  the  main  results  of  this  appendix. 

If  an  exact  linear  relation  exists  between  the  elements  of 

"R-l  / 

V . ±CY  fc, 

h-i  .»j  j k-j 


then 


Some 


m * ° 


In  this  case,  (B-l)  yields 


= X-i  ~ Xk-i-«  = K-i  ' jfi  V'  ^-j  \ 


Therefore  we  can  get  zero  error  by  choosing 


^<r>  J 


0; 

Thus  ^>0  if  an  exact  linear  relation  exists  between  the  elements  of 
*k-i  ) - • ■ > ^fc-p  ■ 

Also  we  have  the  following  general  theorem: 

v-'M 

«=> 


No  exact  linear  relation 
between  elements  of 


HT 


\ 

% K. 


is  positive  definite. 


To  prove  this,  let 


X * 


V 

" J, 
« 

• 

D * 

) 

% 

.V 

a 

. _ 

(B-13) 


(B-14) 


,X$-p) 

( B-l 5 ) 


(B-l 6 ) 


(B-l 7 ) 


(B-18) 


(B-19) 
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Then  F„  s t»M  XK  is  a scalar.  Now  if  and  only  if  an  exact  linear  relation 
exists,  F„  = 0 for  some  D f 0,  no  matter  which  member  function  of  the 
ensemble  we  select  (with  probability  one).  We  also  notice  that 

W - tfXXb  (B— 20) 

and  that  the  ensemble  average  in  (B-20)  is  equal  to  the  matrix  in  (B-18). 

Assume  that  FK  f 0 for  any  D t 0.  Then  /FJ1  >q  for  any  D f 0,  and 
the  right-hand  side  of  (B-20)  is  positive  for  any  D f 0.  Therefore 
XX?  is  positive  definite. 

Conversely  if  is  positive  definite,  the  right-hand  side  of 

(B-20)  is  positive  for  any  D f 0.  Then  |^|*  >0  for  any  D i 0,  yielding 
F*  t 0 for  any  D i 0. 


B-4 
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Appendix  C 


INTERRELATIONSHIPS  OF  Up  AND  Vp 


We  start  with  the  definition  (12)  and  develop  Up  as 


it-S 

.v20C-Wr.-*5X  (by(!5)) 
,vv,^iBTv  (by(12)) 

..^-aW  !bv(13)) 


(by  (13)) 


> Up.,  ~ A^Bp  Up-, 
.(X-  aXH- 


(by  (14)) 


(C-l) 


This  relation  holds  for  p>l,  with  Up=  Rt.  A similar  derivation  for  Vp 


V , p-*i ; V * 


(C-2) 


The  determinant  of  Up  is  given  by 

MV  ub-tfi?)  m 

= MA?  <W«’-ST)  WV  , 


whereas  the  determinant  of  Vp  is 


drV  - *t(r-  8^)  ■ Jet  VM 
- <itl(Ap  -B^)  dft  • <lr{Vp-i 


(C-3) 


(C-4) 


'X  **//*" 


//■  ' ■ 
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Appendix  D 


HERMITIAN  PROPERTY  OF  EXTRAPOLATED  CORRELATIONS 


We  know  that 


We  then  solve 


ft.  = for  lM*f- 


K , i*k*r 


(0-1) 


(D-2) 


for  m , and  set 

1 B / J 


ft?  - f>r  all  for  |k|>  f>-  (D-3) 


We  then  define 


^(p)  _ ^(0  £>r  p4.|  - H • 


(D-4) 


In  a similar  fashion  for  the  backward  case,  we  solve 


- %.  > ) - ^ 1 P 


for  m:  , and  set 


We  then  define 


(D-5) 


ft*  - .4;  W k all  k2i)  R5k,  for  iw*  p-  (D-6) 
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We  know  from  the  definitions  above  that 
Now  we  assume  that 


ff -S?’  W in.  r- 


(D-8) 


^ sKf)  W-  4»ere  f>) 


(D-9) 


that  is,  from  (D-6)  and  (D-3), 


K*  I Ms  I 

Now  from  recursion  definition  (D-6), 

H 


R.(rt  |or  |<k^w. 


(D-iO) 


-«**l  *«  I 

. ±**»  sf 

ft.  I " 

«-i  j=i  nJ 


(%  m) 

(kj  M) 
(lj  (p-i)) 
(ljj  (p-re)) 

(h  (N)) 


(D-1’) 


Therefore  we  have  extended  (D-9)  by  one  step,  and  the  proof  follows  for  all 
kJrp+1  by  induction. 


D-2 
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Appendix  E 


RELATIONSHIP  OF  DETERMINANTS 


The  forward  correlation  recursion  is  given  in  (52)  as 


•d'O  = '4'  a'P’-d  <» 

•'*  “7^*  ’>>»-«  ) - ***• 


The  z-transform  of  this  sequence  is 

ditos  ^2"Rlp)  = i>~A?  . 

M-  | 

The  inner  sum  on  m can  be  expressed  as  (see  (53)) 


OO 


j>M 

Therefore, 


or 


m-±i"  a!?<r>)+ 

&J*  I Vt-  \ 

**■!  ' 


At  the  same  time,  we  define  the  z-transform  of  the  backward  correl 
recursion  as 

k*)  - 

k!*  S PH  * ** 

and  note  that,  via  (62), 


2 2~”  Rlrt  - 

k»*p4-| 


~wt 


(E-l) 


(E-2) 


(E-3) 


(E-4) 


(E-5) 

’tion 

(E-6) 


(E-7) 


E-l 


TR  5501 


A comment  on  notation  is  timely  here.  If  matrix 


£■(*)  5 ?""!)»  , 


where  z is  a complex  scalar  variable,  then 


eH(*) - t'D? . 


(E-8) 


(E-9) 


(E-10) 


which  is  not  always  equal  to  (E-9),  unless  z is  real. 


But  let  us  also  develop  definition  (E-7)  by  means  of  backward  recursion 
(55),  in  a manner  similar  to  that  above  in  (E-l)  through  (E-5).  We  find 

tfv  B f dy  M 


w = i 

The  inner  sum  on  rn  is 


to«p-vi  »*•  i 

-(*-«)  v{p)H\ 


where  we  used  (56),  (2),  (E-3),  and  ( E - 1 1 ) . Therefore 
w'  v ' «*  i r>- 1 


(E-l 1 ) 


VV.  ♦ i = £(?) + f (?) , ,=,« 


/ r i ^ . 

V E-  l J 


(KHW 
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Combining  (E-5)  and  (E-14)  according  to  (E-7),  we  see  that 

d^x-inf)  w(x-  i>~sr)  «• 

must  have  the  same  zeros,  since  these  two  quantities  determine  the  singul- 
arity locations  of  (E-5)  and  (E-14).  The  quantity  $n(z)  defined  in  (E-3) 
is  singular  only  at  z = 0. 

Furthermore 

rip 


±,-Br> , lap.)  , (w, 

where  we  have  utilized  the  observations  that  the  quantities  in  (E-15)  have 
the  same  zeros,  the  same  pole  at  z=0,  and  the  same  scale  factor.  Therefore 
the  two  determinants  in  (E-15)  are  equal. 

Also  since 

det(l  - = 1-2  '4v&,  - - +(-  'f  i^ddSj,  , (e-u 


it  follows  that 


U Af’ . M B f - (cW  By')' 


Numerical  examples  show  that  generally 

|-tr  Af  1 1 \<k*r 


(E-21) 


(E-22) 
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Appendix  F 


SPECTRUM  FROM  EXTRAPOLATED  CORRELATIONS 


The  forward  correlation  recursion  is  given  by 

r'  > 


where 


r 

i 

i 


, hii  ? 

and 

We  wish  to  evaluate  the  z-transform  of  •' 


oo 


>lp) 


A*®-'**  W 

In  order  to  do  so,  consider  a fictitious  process  [X„] 
lation  given  by  (F-l;  through  'F-3  . Consider  the  output 
predictive  error  filter,  jiven  py 


A? >L , “I'  h 


The  crosscorrelation 

c. . %U.  • - ■ - iATR:: 

’ *.o  *••• 

using  7 and  **  *',a* 

A A 

(^  * 0 fcr  * j (C.fC 


JE= 


,1 


(F-l) 

(F-2) 

(F-3) 

:F-4) 

with  the  corre- 
of  the  optimum 


I » . r 

, all  » 

o), 


V 
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that  is,  predictive  error  filter  output  Y„  is  uncorrelated  with  all  past 

A 

values  of  input  )(*  . 

Also,  output  autocorrelation 

l-TZ-Z&Z  Ar“-£e~£\ 


(F-8) 


using  (F-5)  and  (F-6).  But  now  employment  of  (F-7)  in  (F-8)  shows  that 

3)^=0  (F_9) 

Also  (F-7),  (F-6),  (F-2) , and  (12)  yield 
And  since,  from  definition  (F-8), 


(F-10) 


a H 


X-XL  , 


(F-ll) 


we  have 


Uf,w  = o ") 

0,  o&enHjtj  * 

A A 

that  is,  predictive  error  filter  output  ^ is  white  for  input  X(<.  (Of 
course.  Up  is  not  diagonal). 

A 

At  the  same  time,  autocorrelation  Dm  can  be  expressed  (by  means  of 
(F-5))  as 


(F-12) 


A.  ■ £ K < - % | *?<■  Af  . •»  - 


(F-13) 


Therefore  the  z- transform  of  {*  is 


F-2 
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±rt>.  • * £****0*.  ±2' a;,,h 

t»».<w  »l=0  !»«•<«  js-0  J 


where  we  have  used  (F-13),  (23),  (F-4),  and  (E-9).  When  we  couple  (F-14) 


(F-14) 


with  ( F-l 2) , we  obtain 


- % to  Jf  QlZ'V) 


( F-l 5 ) 


( F-l 6 ) 


where  matrix  Up  is  independent  of  z.  This  is  one  of  the  main  results  of 
this  appendix. 


If  we  let  (for  f real ) 


•=  exp(i?irfa)  ; )f) 


(F-17 ) 


< 24  ; 


and  denote  the  forward  predictive  error  filter  transfer  function  and  spectral 


estimate  as  - - ±*& ■ hTw. 

respectively,  then  the  spectrum  of  process  can  be  expressed  as 

G?’(f)  - A W V)p  , 

where  we  have  utilized  the  result  that  (see  ( E-8) ) through  ( E- 10) ) 


( F-l 8 ) 


(F-19) 


**"*  *-1**$t * £ *»  %+*  -<**.**  ^Si  ** 
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The  procedure  for  the  backward  correlation  recursion  parallels 
above  to  yield  (using  (23)) 


Opu)  - 

*C7 


w - * Kwf  v,  fj«" •(,-)] 


and 


-i 


<(frvP  Hr^)"H . 


r 


ZERO  LOCATIONS  OF 

Assume  that  ^(3)  1 c Q(z)  has  a zero  at  z=zf  ? 0;  that  is, 

ftssu*e  Q(?|)  = °)  "tW  0 

where  0 < |?,|  . But 

nfa  - - if  /d = - if[^  *c  + • + rtf-  / r 


Therefore  iff  fa)  is  finite  for  0<  |e,|  , yielding 

Q U)  ?£’(?,)  = o * i 

Therefore  assumption  (F-23)  is  invalid,  indicating  that 


Q(S-)  +0  for  0 ^ 1*1  . 


Now  from  (F-24) 


-*rtf  « N-»  Pj 


F -4 


that 


(F-21) 


(F-22) 


( F-23 ) 


(F-24) 


(F-251 


[F -2b  ) 


’*-2" 
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therefore. 


QU)-  - 2f  ^ to  0< 


( F-28 ) 


Thus  Q (z ) has  a p-th  order  zero  at  z=Q,  but  is  not  equal  to  the  zero  matrix 
for  0<N  . Of  course,  the  individual  elements  of  matrix  Q(z)  can  have  zeros 
anywhere . 

POLE  LOCATIONS  OF  ' 

n 

Since  from  (F-24) 


(F-30) 


(F-29) 

where  Qp(z)  is  a matrix  of  polynomials  in  z of  order  p,  it  follows  that 

QW=-*fOP(*)"=-  (f-30) 

A 

where  <Ww  is  a matrix  of  polynomials  in  z of  order  (M-l)p.  Therefore 
the  poles  of  Q(z)  are  caused  by  the  zeros  of  det  Qp(z);  that  is,  the  poles 
of  W are  caused  by  the  zeros  of  det  . As 

therefore,  Q(z)  ~ I as  so  that  Q(z)  has  no  poles  at  !*)  = «• 

Thus  the  poles  of  Q(z)  are  located  where  det  ^V)  - o. 

We  now  consider  the  problem  of  determining  when  det^]^(»)rO:  the 


following  derivation  is  based  upon  Ref.  7.  Let 

,/>■] 


be  an  Mpxl  matrix.  Define  prediction 


(F-31) 


i - c 1, , 


( F-32) 


F-b 
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where  C is  MpxMp.  And  define  error 


Then 


where 


K-t-i-X-ci.,. 

uf  • {X-  C !)(?"- 1!  c*5 
* % - c*<_- ft,  c"+  c%  cH 


(F-33) 


(F-34) 


-\H 


(F-35) 


u 

The  minimum  value  of  tr  is  realized  when  (see  appendix  B)  we  select 


uk 


(F-36) 


The  corresponding  value  of 

S,  C = % -%KX  - % - cyv  cH , (|;-37> 

Ij 

since  Now  let  the  left  eigenvectors  and  eigenvalues  of  the  optimum 

C be  denoted  as 


?”C  = V.?*  > 'iWS  MP 


(F-38) 


(The  eigenvectors  may  not  all  be  linearly  independent).  Then 


0^?:  s„f  = £ s.  £ t -£{iL-ai  cH)l 


(F-39) 


F-6 


Now  % is  Hermitian,  block  Toeplitz,  non-negative  definite,  and  has 

X •••  v 

V-- 

Kf  V 

Therefore  |\J<|  for  that  is,  all  the  eigenvalues  of  C are 

bounded  by  unity  in  magnitude.  Furthermore,  Ref.  7,  p.  134,  shows  that 
if  there  is  no  exact  linear  relation  between  the  elements  of  X„  X*.,  ... 
then  j V^J < | for  |<w*Mp  (see  also  appendix  B). 

Now  we  develop  the  error  in  (F-33)  in  more  detail: 


Wl  l\.  I Wl  III 


V.  - II 


r 

>s  ■ 

1 

c.  •• 
« 

■ c; 

r 

T 

LVJ 

■V 

J 

I/'** 


Mi 


inimizing  tr  can  be  seen  to  make  C of  the  form 


C - 


aT  Af  a; 

x o o 

o r o 

. % 

\ . 

O TO 


V' 


Therefore  (Ref.  7,  eqs.  (35)  and  (36)), 


da(c-xr)  = (-  xf P ^\\) 
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If  we  were  to  assume  that  where  l*t|£  I , we  would  have  det(C-*,l)~  0. 

But  this  contradicts /X*}<  I for|<*sMp.  Therefore,  the  zeros  of  det  %%) 
all  lie  inside  the  unit  circle;  that  is,  the  poles  of  Q(z)  all  lie  inside 
the  unit  circle. 


F-8 
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Appendix  G 

HERMITIAN  PROPERTY  OF  ONE-STEP  EXTRAPOLATED 
CORRELATION  MATRIX  ESTIMATES 


From  (78),  at  the  (p-l)th  stage,  we  know  that 


SlC  •S'lC.Bt’"-  f aLV.  * K , u«  s f-1- 

M«I  «-«t  •»  «.  i mi  7 


(6-1) 


Now  we  start  with  (94)  and  express 


SlxIX  - H#'" 

M:l  T » r 


(h  (M) 

fa  w) 


(6-2) 


Thus,  the  one-step  extrapolated  correlation  matrix  estimates,  based  on 
order  p-I,  are  Herrmtians  of  each  other. 


G-l/G-2 
Reverse  Blank 
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Appendix  H 


INTERRELATIONSHIPS  OF  U AND  V FOR  UNKNOWN  CORRELATION  CASE 

P p 


We  develop  the  definition  (S5)  as  follows: 


Ir  jv 

i 


’©4 


I 

>1 


> * 
f 


1 €>, 


Now 


Therefore 


4 = 

.-2./tX*$£i£iL. 

M*<>  ' 

(M™)) 


%$X.r  - 

r»«  1 r j*.i  ' J J 

= vr,  .‘ftSV 


'-P 


(i>j  ()*)) 


(j3bs)) 


In  a similar  manner,  we  can  show  that 

V(:n-B?AT)Vr, 


(H-l) 


(H-2) 


(H-3) 


(H-4) 


In  order  to  show  that  Up  is  Hermitian,  we  recall  the  constraint  (98) 


i: 

i! 


\ 

1' 


u 


s 

i 

1. 


I 

< ■■v. 

!<• 


and  express 


$ 


H-l 
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i| 


H 


*.  * 
V>v  f 

•] 


* 


% 


0rVJr-A^Ur-Ur-.fV;,A7. 


#>>„>■  .if)" 


(H-5 ) 


Therefore  if  Up-i  - ULj  *'  . it  irmiediately  follows  that 


up  . uP 


( H— 6 ) 


Similarly  since 


it  also  immediately  follows  that 

r« 


v;-vr 


(H-7 ) 


3ut  properties  (H-6)  and  (H-7)  are  obviously  true  for  p = 0,  because 


(H-8) 


J0  *6  *0  o 

Therefore  \H-6)  and  (H-7)  are  true  for  all  p,  by  induction. 


In  order  to  relate  the  determinants  of  Up  and  V , we  express  (H-3) 


and  (H-4)  as 


Ur  < (Af- B^,  , V,  ■ 


(H-9) 


Therefore  if  det  Up- 1 = de-fc  Vjs-/) 

Up  = drt  Vp  • (H-10) 

But  (H-10) is  obviously  true  for  p=0  by  (H-8).  Therefore  (H-10)  is  true 
for  all  p,  by  induction. 

Prope-ites  (H-6),  (H-7),  and  (H-10)  applied  to  (98)  immediately  show 

that 


(H-ll) 


H-2 


i; 

?■ 


V 

1 ■ 


1i 


k 

ft 

i 


ii 


II 


A . 
? - 


5 

r 


* 

? - 


j; , 

'i  • 
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.*3 
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Appendix  I 


MINIMIZATION  OF  TRACE  OF  WEIGHTED  ERROR  MATRICES 


We  wish  to  minimize  the  trace  of  (123)  by  choice  of  matrix  Gp.  We 
use  the  fact  that,  for  square  matrices  P and  Q, 


4r(PO)  = 2-P.C..  20- P„-  tr(QF), 


(1-1) 


to  express 


{r  Ep + rr  r “MX-.  S'r  + 5^  - <V^-.  ^ fj ^ V'  '4-' 


(1-2) 


Now  (1-2)  is  an  analytic  function  of  the  variaoles  Re(Gffln)  and  Im(Gmn). 
Therefore  the  minimum  of  (1-2)  is  realized  simply  by  setting  the  coefficient 
of  Gp  equal  to  zero  (Ref.  20).  We  obtain,  after  premultiplying  by  and 

post-multiplying  by  , the  equation  for  G^: 


c+a"  u;,1  ■ s;tvpv,  £ +jv,  ip:  sp  (1.3) 


(Gp  is  not  Hermitian  or  Toeplitz,  as  numerical  examples  will  show.) 
terms  of  A^and  B^5,  we  have  the  simultaneous  equations 


where  we  utilized  (122). 


In 


(1-4) 


I-1/I-2 
Reverse  Blank 
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Appendix  J 


COMPUTATION  OF  FILTER  TRANSFER  FUNCTION 


The  forward  predictive  error  filter  transfer  function  is  given  in 


(68)  as 


Hf(0  = - ^ exp(-;  t |f| 


divide  the  frequency  range  (^j^-J  into  Np  cells  of  width 


A = J-  -L 

^ A 


Then  for  |»»t|  5 


Ha  (w  - Ha  (^4)  = - £ ?*?(-'  A1? 

|k 


(J-l) 


(J-2) 


(J-3) 


£ 0;  pf-IS  " S 

Now  if  we  let  the  sum  in  (J-3)  be  denoted  as  an  Np-point  FFT, 


S*L  0 s * s H-1/ 


then  (J-3)  becomes 


*(*■)■ 


2*,  ; 0SN<  NV/2 

, -hJ2  i ws-i 


Then  quantity  ^ in  (J-5)  is  an  MxM  matrix  for  each  value  of  m. 


(J-4) 


(J-5) 


(J-6) 


J-l /J-2 
Reverse  Blank 
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Appendix  K 


PROGRAM  FOR  SPECTRAL  ANALYSIS 


In  this  appendix  we  present  the  program  for  tne  procedure  summarized 
in  (184)  - (187)  and  (165).  The  spectral  estimate,  (165),  is  computed 
at  frequencies  fa/(Nf6)]  ' 


[)r  hT(^9  > *"•» s h/2; 


(K-i) 


where  the  forward  predictive  error  transfer  function  is  given 

by  (J-6).  The  specific  scalirg  adopted  is  based  upon  (166),  which  takes 
the  sampled  form 


(K-2) 


where  [wm}  is  a set  of  integration  weights  (e.g.,  trapezoidal).  The 
approximation  is  a good  one  if  G^(f)  is  sampled  finely  enough;  that  is, 
if  Np  is  large  enough  to  resolve  the  peaks  and  valleys  of  G^P)(f).  If  we 
employ  (J-2),  (K-2)  becomes 


( K-3A) 


or,  for  trapezoidal  weighting, 


(K-3B) 


where  we  have  employed 


the  periodic  nature  of  G^(f)  (See  (165)  and  (68)). 


K-l 
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Thus  the  sura  of  samples  equals  the  sample  power,  (80). 


For  a real  multivariate  process,  we  can  employ  (171);  a modified  form 
emerges: 

Utfi 

^ frQ  *-  5 ^ ^ r ^ p*0*"'  (k-4) 

where  |55^  is  another  set  of  integration  weights.  This  is  the  form 
programmed  in  the  following;  the  quantities  computed  are 


The  real  part  of  their  weighted  sum  equals  the  sample  power,  Rq.  The  FFT 
used  here  to  evaluate  (J-5)  is  given  in  Ref.  21;  it  is  limited  to  powers  of 
2,  hut  could  be  replaced  if  desired.  Input  parameters  are  N,  PMAX,  and 
NF  in  line  22,  and  the  input  data  call  is  in  line  37  and  SUBROUTINE  DATA; 
all  these  quantities  have  to  be  changed  by  the  user  to  fit  his  particular 
application.  The  program  is  written  for  a real  multivariate  process 
(general  M),  with  the  exception  of  FUNCTION  DETERM,  SUBROUTINES  SDM, 
INVERT,  and  SOLVE,  and  the  printout  of  the  spectral  density  matrix,  (K-5). 
Arrays  used  in  the  program  are  explained  by  comment  statements.  A sample 
printout  follows  the  program. 


K-2 
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c MuLTIYAHIATt  LINEAR  PREDICTIVE  SPECTRAL  ANALYSIS* 

c Employing  stiGhTtu  forward  and  backward  averaging, 

C This  PROGRAM  is  MITTEN  for  RE*L  PROCESSES  and  general  M»  *ITH  Thr: 
c Exception  of  function  determ  and  subroutines  som.  invert*  anc  sol/E, 
C Auj  I HE  print  cut  of  the  SPECTRAL  DENSITY  matrk. 

C USEh:  change  LINES  22  AND  37*  AND  REPLACE  SUBROUTINE  OATA, 

C m r jlMtNblONALlT Y OF  MULT I VAR J ATE  PROCESS*  INTEGER  INPUT 
C N = NOKdth  OF  DATA  POINTS  IN  EACH  PROCESS*  INTEGER  INPUT 
L X ( 1 # i ) , , , X ( N * 1 ) » • . , * X ( 1 * M ) , , . X ( N ? M ) = INPUT  LATA*  ALTERED  ON  OUTp-iT 
C PhAx  = MAXIMUM  ORDER  OF  FILTER » INTEGER  INPUT 

l  NF  = SIZE  OF  FFT  (MUST  dE  A POnER  OF  2 TO  USE  M*LFFT)j  INTEGER  INrUT 

C AVE  z MEANS  OF  INPUT  uATAl  OUTpUT 

C R : COVARIANCE  MATRIX  OF  INPUT  DATA  I OyTPUT 

C AlC  z AKAIKE’S  INFORMATION  CRlTERIONl  OUTPUT 

C PbEsT  = BEST  ORDER  OF  FILTER*  INTEGER  OUTPUT 

l UoESI  = MATRIX  OF  COEFFICIENTS  IN  SPECTRAL  ESTIMATE?  OUTPUT 

t AP  = MATRIX  OF  FOK*»ARU  PARTIAL  CORRELATION  COEFFICIENTS*  THEN  = 

L MATRIX  OF  FORWARD  PREDICTIVE  FILTER  COEFFICIENTS  FOR  PBEST*  OUTPUT 
C dp  = MATRIX  OF  BACKWARD  PARTIAL  CORRELATION  COEFFICIENTS*  OUTPUT 
l XX*YY  = SPECTRAL  MATRICES*  OUTpUT 

PARAMETER  M=2  U BIVARIATE  PROCESS 

PARAMETER  N=  100  » PMAX=  10*  NF=102r»  nF41=nF/4+1 

INTEGER  PdEST »P 

DIMENSION  X(N»M)  »Y(N»M)  ,Z(N,M>  *UBEST(M*M)  , Ap  (N»M,p;‘,AX ) * 
4LP(M*M,PMAX) iAVE(M) *XX(NF»M,M) *YY(NF»M*M) *CCSI (NF41) * 

$U  (M»M)  *V(M*M)  »UI(M»M1»VI(M*m)  » A { M » M ) »B(M*N)  ( r ( ,-1  * M ) t 
Sv.A(M*M)  ,Wd(M»M)  »WC(M»M)  *WD(M»M)  »WE(M»M)  *AIC(P*>  AX)  , AIC0(2) 
EQUIVALENCE  (X*Y) * ( AIC ( 1 > »AiC0(?) ) 

C PRINT  OUT  VALUES  OF  PARAMETERS 
I=N 

J=PMAX 

K=M 

L-NF 

PRINT  1,  I » J * K * L 

1 FORMAT ( 1H1 * * N =* * I6*10X* »PmAX  r**l4»19X*’M  = ' r 12  * 10X * *NF  =*»I5) 

C INPUT  DATA  IN  X ( 1 * 1 )*.. X (N» 1 )*,*.# X ( 1 »M) ... X ( N, M ) 

CALL  DATA 
PRINT  2 

2 FORMAT (/'  INPUT  DATA:*) 

J=N-99 

L=N-200 
DO  3 I=1*M 
PRINT  4*  I 

IFUM.LE.200)  GO  TO  5 
PRINT  b*  (X(Kr I) >K=1#100) 

PRINT  It  L 

V FORMAT ( lb*  * INPUT  DATA  POINTS  NOT  PRINTED  HERE’) 

PRINT  o*  (X(K* I) *K=J»N) 

00  TO  3 

S FRINT  t>t  (X(K*I)  *K=1*N) 

3 CONTINUE 

4 FORMAT ( * PROCESS  NUMBER* *l2) 

b FORMAT (SE20.8) 

C EVALUATE  PARTIAL  CORRELATION  COEFFICIENTS 
CALL  PCC 
PRINT  8 


K-3 
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o f OKMAT ! / ’ MEANS  OF  INPUT  DATA!’) 

PRINT  b,  (AVE(I) *I=i»M) 

PRINT  9 

9 F OKMAT (/•  COVARIANCE  MATRIX  OF  INPUT  DATA:*) 

PkAimT  o*  ((R(I»,J)»I=1*M)»J=1#M) 

PRINT  10 

10  P OKMAT ( / * AKA1KE  INFORMATION  CRITERION:*) 

PRINT  11»  (P#AIC(P) »P=0,PMAx> 

11  FORMAT ( I10»E20«d) 

PRINT  12»  PBEST 

12  FORMAT!/*  PBEST  =*»I3) 

PkINT  13 

13  *•  OkMAT  (/•  UbEST:*) 

PRINT  b,  ( (UBESTd*  J)  , 1=1, M)  , J=1#M1 
PRINT  14 

.14  f OKMAT (/•  FORWARD  PARTIAL  CORRELATION  COEFFICIENTS:’) 

L 0 lb  P=1*PMAX 

lb  PRINT  lb  , P,  ( !AP(I,  J,P) , I = 1,M)  , J=ltM) 
lb  F URMAT(I10»bE20.a) 

) HINT  17 

17  FOkmaT!/?  bACKrtARD  PARTIAL  CORRELATION  COEFFICIENTS  J ' ) 

LO  id  P=i*PMAX 

18  PRINT  lb*  P»  ( (bp  ( I * J#  P)  , I“1  f M)  * .J=l  * M) 

IP (PbEbT^LO.O)  GO  TO  19 

L EVALUATE  PREDICTIVE  FILTER  COEFFICIENTS 
CALL  PFC 
PRINT  20 

20  FORMAT!/’  FORWARD  PREDICTIVE  FILTER  COEFFICIENTS  FDR  PbEST : * ) 

LO  21  Prl*PBEbT 

cl  PtUNT  16*  P, ((AP(I»J,P)fI=l,M)»J=l*K) 

C EVALUATE  PREDICTIVE-EkROR  FILTeR  TRANSFER  FUNCTION 

19  CALL  PEFTF 

c evaluate  spectral  density  matrix  and  coherence 
K=NF/2+l 
call  sum 
PRINT  22 

€.2  FORMAT!/*  SPECTRAL  DENSITY  MATRIX  AND  COHERENCE  FuR  M=2:») 

PRINT  23 

c3  FORMAT  (8X*  ’BIN*  *10X,  JAUT011**14X*  *AUT022’  *.10X,  » HEAL  (CROSS12)  ’ , 7X  * * 
$1iv.aG(CKC’SS12)  ’ *9X* ’MAG  SQ  CqH’  * 11X*  ’ARGUMENT  * ) 

PRINT  16*  (L»XX(L»1,1),XX(L,2*2) »XX(L,1,2) ,YY(L*1*2) »YY(L,1,1) fYY( 
$L»2*2) * L=1*K) 

SUBROUTINE  DATA 

C THIS  SUBROUTINE  GENERATES  DATA  FOR  M=g»  BIVARIATE  PROCESS 

DEFINE  IRAND?1*5**15+ < ( 1-SlGN ( 1 * 1*5**151 ) /2 ) *343597383b7 
DtF INE  RANp=F;LOAT  ( I ) /3H359738367 , 

1=5281 

TA=0. 

TB=0, 

DO  1 KsltlOO  iD  WILL  DISCARD  THESE  INITIAL  POINTS 
I=IRAND 

Tn.85*TA-.75fTB+RAND»,5 

I=IRAND 
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To= . 65*T A* . 55*Ta+RAN0- , 5 

1 TA=T 
Xil#l)=TA 
aU»2)=TB 
DU  2 K=2#N 
UlRAND 

T=  . 85*T A- . 75»T0*RANO*  . 5 
1-1RANO 

To= . 65*T A* . 55*  TBtRAND- , 5 
Ta=T 

XIK,1)=TA 

2 XlK»2i=TB 
hETURN 

SUUhOUTINE  PCC 

L T.ili,  SUrtkUwTlNE  COMPUTES  PbEST,  UbEST*  AND  THE  PARTIAL 
t CoKNELATION  COEFFICIENTS  FOR  P = i TO  pMAX  f ANY  F 
1-n 

vi=PMAX 

lA=0.*b0RT(N)/M 

lr  IPMAX.GT.IA)  PRIM  1,  J#I,IA 

1 FOKMAT(/'  PMAX  =*#I4,»  IS  TOO  L«RGE  FOk  NUMBER  OF  DATa  POINTS  fcl  =* 

SEARCH  LIMITED  TO  P =»»14) 

1 A— M I N { I A » PMAX ) u)  UPPER  BOUND  ON  PMAXJ  E;  163 

F«L=2.***M/N  'J  FAr.=0,  WOULD  FORCE  PB^ST  t JliAL  TO  PMAX 

c suotkact  means * fill  in  data  arrays;  eq  no 

LO  2 I - 1 » M 
TA=0, 

uU  3 K=1#N 

3 I A=TA+Y(K» I ) 

ta=ta/n 

AvE(I)=TA 
UU  2 K=1»N 
Y(rwI)=Y(Kf 1)-TA 
C 2(K»I)=Y(K»I) 

L INITIALIZE  CORRELATION  MATRICES*  EOS  62 1 114,  Ak'L  10b 
CmLL  AUT0(2»N-1»  Y,  ,iC) 

10  4 1-1  > M 
DO  4 -J=I»M 
TA-Y(1, I)*Y(1» J) 

Tl>=Y(N,I)*Y(N»J) 

K(I»J)=(WC(I,0)+TA+TB)/N 
w A ( I » J)=WC(I,  J'y+TB 
V>U(I»J)=WC(I»J)+TA 
K{vl»I)=R<I»J) 
wa(J»I)=wA(I,J) 

4 AblJ»I)=Wb(I.J) 

CALL  EQUAL{R»U) 

CALL  EUUAL(R»V) 

CALL  CKOSb(2,N,Y»Y»WC) 

C BEGIN  RECURSION 

AIC(0)=LOG(DETERM(U) ) 

AiLMINzAlUO) 

FBEST=0 
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CALL  EQUAL (U#UHESTJ 
Do  5 P=1»IA 

L EVALUATE  MATRICES  REQUIRED  IN  BILINEAR  MATRIX  EQUATION!  tQ  126 
CALL  lNV£RT(V#VI) 

CALL  MULT (VI » kd» WD) 

CALL  EuUAL(*D#«fd) 

LmLL  InvERT(U»UI) 

CALL  EQUAL (mA**D) 

CALL  MULT(WD#UI#ttA) 

CALL  AUDC*C'0C#mC) 

L SolVE  BILINEAR  MATRIX  EQUATION;  EQS  157-161 
CmLL  solve 

I evaluate  partial  correlation  coefficients!  eq  i^c 

call  MULT(i*C#VI#A) 

CALL  TRANSUC*»D) 

CALL  MULT(WUrUI,B) 

CALL  EuUAL(A»AP(l,i,P) ) 

call  Equal (b»3p ( l » l » p ) ) 

C UPUATE  MATRICES  U AND  V ! EQ  18i 
CALL  MULT(A»WD»V*E) 

call  Sub ( Uf we* u) 
call  mulT(b»mc^E) 

CALL  SUB  CV»hE» V) 

C CALCULATE  AKAIKE'S  INFORMATION  CRITERION!  EQ  180 
A I C ( P ) =LOU { DETERS ( U ) ) +F AC*P 
lMAIC(PJ.GE.AICMIh)  GO  TO  6 
AICMIN=AIC(P) 

PdEST=P 

CALL  EQUAL{U#UBEST) 
b IF (P.EQ, IAJ  GO  TO  5 

C UPDATE  DATA  SEQUENCES  Y AND  Z 1 EQ  111 
L=P+1 

LO  7 K=N»L»-1 
Lu  8 I — 1 » M 
TA=Z(K-1»I) 

CO  9 J=1#M 

9 TA=TA-B(lrJ)*Y(K> J) 

b c(K»I)=TA 

DO  10  I=1>M 
T A-Y (K » I ) 

LU  11  J=1#M 

II  T«=TA-A(I»J)*Z(K-1»J) 

10  Y (K » I ) =TA 

7 CONTINUE 

c calculate  new  correlation  matrices!  eq  n+ 

CALL  AUT0(P+2>Ni Y»*A) 

CALL  AUTO(P+l>N-lfZ,WB) 

CALL  “ROSS ( P+2 »N»YfZfWC) 

5 CONTINUE 

IF(M,EQ,1)  RETURN 
K=M-1 

DO  12  1=1 t K 
L=i  + 1 

DO  1?  *J=Lfrt 


¥ 
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DdtbT  1 1 • J)=#5* (UBEbT  ( I » J)  ■MJhEST  ( J»  I)  ) 

12  ubLbr<j,i)=uBLbni#j) 
return 

SJuKOUTINk.  PFL 

C Tdlb  bUbKOuT iNE  CG.iPUT  Eb  TrtE  PREDICTIVE 
L ULTcK  COEFFICIENT ANY  m;  td  79 
lFCPBEbT.LE.l)  RETURN 
CO  1 P=2  »PbEST 
1A-.P-1 
LO  2 L-l * IA 
Io=P-L 

caLL  MulT(AP(1»1»P)  »bP{  i , 1 , jB)  » <iA) 

CALL  bUd(AP(l»l,L)»VkA,WA) 

Call  MULribP(l#l»P) ,AP(i,l,L) ,fcU) 

CALL  SUB  (lip  ( 1 » 1 » lb)  #WB»bP(l,l»Ih)  ) 

£ CAlL  tu«  IAL  ( WA » AP  ( 1 » 1 » L)  ) 

1 CONTINUE 
Kl 1 URN 

L 

bUoROUIJNE  PEFTF 

C T.tlb  SUBROUTINE  COMPUTES  TuE  PrEDICTI Ve-EKHOK 
^ PIlTlH  TKhIISFEK  FumCTION)  ANY  Ml  eOS  bB  AND  iJ-i)-(J-o) 
l\::l.L427*L0G(NF)+.b 
CALL  dTRCOS(COSl»NF  ) 

L O 1 I - 1 » M 
Lu  1 J=1*M 

xx(i»i»j)=o. 

IMI.E«I,J)  XX(1«1»J)=1. 

Yr  (l»l»vJ)=0. 

IF  (PbEbT/EQ.O)  t>0  TO  2 
l«=PbEST+l 
DO  3 L=?»IA 
XX(L*I*J)=-AP(I,J,e-1) 

0 YY(L*I»J)=0. 

2 Im=PUEST+2 
lu  A L=IA»NF 
XX(L» I t J/=0, 

9 YY(Lil»J)=0, 

1 Call  MKLFFT(XX(i»1»J)»yy(1»i»J),COS»I»K»-1) 

Re  TURN 

c 

bUbHOUTINE  SDM 

C Thlb  bUBKOl iTINE  COMPUTES  TnE  SpECTRAL  DENSITY 
C MATRIX  AND  COHERENCE  FOR  M=2J  EQS  178  AND  K-L 
T=2./NF 
DO  1 L-l » K 
ViA  ( 1 r 1 ) =XX  (L  r 2 r 2) 
wA(l»2)=-XX(L*l»2) 
i\A  (2»  1 ) e-XX  (L*  2*  1 ) 

*.A(2»2)=XX(L»1»1) 
i\o C 1 # 1 ) =YY  (L»2»2) 
hd(l>2)E-YY(L»l»2) 

;.d(2*l)=-YY(L»2»l) 
aB(2»2)  = TY ( L » 1 » 1 ) 
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T *=0ETERM  { WA ) -OE.TE.KM  ( WB ) 

Tb=t»A  (lfi)«MB(2»2)twA(2>2)*wB(l»l  )"WA  (1#2)*alm2»1) 
TA=T/(TA**2*Tb**2) 

CALL  tkansua#»c) 

CALL  MUlT<UBESTf*CfWU) 

CALL  MULTUB»WJ»HC) 

1 o=rtC  ( 1 • 2) -WC  (2»  1 ) 
call  MUlJ  (WAf*UfWC) 

CALL  TRANSl»B»wO) 

Call  mulTI UbESTFWDFwE) 

Call  mult  UbfWEfWOJ 
CALL  A00{*»C#WU»«C) 

YY  (L#1»1)  = (WCUf2)**2+Tb**2)/(WC(1#1>**C(2»2)  ) 
f t IL*2*2)=A1AN2(  Tbp»Ca,2)  ) 

*X(Lf1f1)=TA*wC(1»1) 

aX  l L * 2 f 2 ) -TA^fcC  (2 f ) <*• 

XxILf1f2)-TA*WC(1fl) 

Y Y(L»1»2)=TA*TB 
XX(L»2# 1 )-0. 

Y Y IL » 2 » 1 )-0» 

1 CONTINUE 
hL l URN 

C 

SUBROUTINE  CR0SS(N1.N2#A#B»c)  0 A>BfA  NG 
C Tulb  SUBROUTINE  COMPUTES  A CROSS  CORRELATION  PATH IX! 

L IMENSICN  A(N»M) fB(N»M) ,C(M,M) 

DOUBLE  PRECISION  T 
L-0  1 1 = 1 fM 

LO  1 J=1 f M 
TrO.UO 

LO  2 KSN1’N2 

2 T=T+A(K»I)*b(K-l»J) 

1 C(I»JJ=T 

return 

c 

SUbROUTlNE  AUTO ( N1 f N2 f A f B ) » AfA  NG 

C Trilb  SUBROUTINE  COMPUTES  AN  AUjO  CORRELATION  MATRIX} 
DIMENSION  A(N#M) Fb(MfM) 

LOUdLE  PRECISION  T 
LO  1 I = 1»M 
LO  1 J=I#M 
T=0.DP 

LO  2 K=N1fN2 

2 T=T+A(KfI)*A(KfJ) 

1-  ( 1 f J ) =T 

1 h<JfI)=B(IfU) 

rlTurn 
C 

SUBROUTINE  EQUAL(AfE) 

C This  SUBROUTINE  SETS  TWO  MXM  MATRICES  EQUAL 
DIMENSION  A(MfM) Fb{M»M) 

LO  1 1 = 1 * M 
CO  1 J=1»M 
1 R J)=A (If J) 

RE  I URN 
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-fcA(2*l)*«B, |f?) 


MAo  SC  COH 
argument 
AUTOll 
AUT022 

REAL (CROSS] _ ) 
IMAG( CROSS 1 ) 


ANY  4l  LQ  11 *R 


ANY  Ml  EG  11., A 
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J 


■ * 

4 


SUBROUTINE  TRANS(A*B)  fi  A, A NG 

C Trtlb  SUBROUTINE  TRANSPOSES  AN  MXM  MATRIX 
t-iMENSlON  A(M#M)  »B(M#M) 

LO  1 I=1#M 
uO  1 0=1 »M 
1 t tI*J)=ACJ»l) 

F E T URN 
c 

bUbROUTINE  ADU(A*b*Cl  ft  A#ri*A  OK 

c inis  subroutine  auos  two  mxm  matrices 

LlMt.'JSlCN  A(M#M)  ,C(M,M) 

LO  1 1 = 1 » M 
LO  1 J=1»M 

1 C (1)  J)=A(1»J)+B(1»U) 

RETURN 

c 

SUuROUTINE  SUb(A#b»C)  it  A »n*  A OK 

C Thlb  SUBROUTINE  SUBTRACTS  TWO  mXM  MATRICES 
LiNENSlON  A(M»M) #b(M,M) ,CIM,M> 

L’U  1 1 = 1 * M 

CO  1 J=1  »M 

1 C(i#J)=A(IrJ)-U(I#J) 

KtTURN 

C 

bUbROUTINE  MULT(A»d*C)  U A,B»A  NG 

C Thlb  SUBROUTINE  MULTIPLIES  TWO  MXM  MATRICES 
L i MENS I ON  A (M»M) ,b(M#M) ,C(M,M) 

LO  1 1=1 »M 
LU  1 J=1»M 
7=0, 

CO  2 K=  1 » M 

2 T=T+A(I,K)*b(K» J) 

I C ( A » J) =T 

RETURN 

C 

bUbROUTINE  INVERT  { A>B)  fu  A t A NO 
L Thlb  SUBROUTINE  INVERTS  A 2X2  yATRIX 
DIMENSION  A(2f2) *B(2>2) 

TA=1./UETERM(A) 
i>(1,1)=A{2»2)*TA 
b(2,2)=A(l#l)*TA 
b ( 1 » 2 ) =-A ( 1 # 23  * TA 
b(2, 1)=-A(2#1)*TA 
KETURN 
L 

bUbROUTINE  SOLVE 

C THIS  SUBROUTINE  SOLVES  blLlNEAR  MATRIX  EQUATION 
L FOR  -|=2»  BIVARIATE  PROCESS?  EQ$  1S7#  1&8,  ANB  182 
TM=WAa,l)+»vA(2»2)+WB(l,l)+wB(2»2) 

TB=UtTERM l W A ) -UETERM ( Wb ) 

CALL  MULT(WC»WB»WD) 

«vt_  ( 1 * 1 ) =rtA  (2»2) 
ac.  ( 1 * 2 } =-WA  ( 1 » 2) 
ac(2»1)=-WA(2»I) 

At ( 2 * 2 ) rwA  del) 


I 

1 

l 

l 


1 

i 

i 
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CmlL  MUl.T(V»E»fcC»!«A) 

CALL  ADL'U&#WA.kO) 

Add#  1 ):TA*^B(i»l)tTB 
*D(2»2}~TA#ftBI2#2)+Tb 
*cHi*2)=TA**Bll*2> 

,.ii(2#l)=TA»WB(2tl) 

Call  INVERT { Wb* wEi 

call  MULT(WD»WE»WC) 

KtTUKN 

c 

FUNCTION  DETEKi«1(A) 

c this  (-unction  computes  the  determinant  of  a 2x2  natrix 

L ILLUSION  A(2»2) 

Cl  I lHMsA {l»I)*A(2»2)“A(l»2)*A(2»i) 

i-t_T<JKN 

t.  IMU 

EouHOuTINE  MKLFFT(X*Y,CC»M»iSN) 
l IMENSION  X ( 1 ) »Y(1) ,CC{1) »L(12) 

equivalence  (Li2»ui>>»  (H1,l<2)  >»(Uo*l(3)  ) » il9»l(4)  > , (lb»l(5)  ) # 
1 (*.7»L(o)  ) » (L6»L(7)  ) » (L5»L(8)  ) » (L4#L(9)  ) # (L3»L(10)  ) » (L2»L(11)  ) » 
2(cl»L(12) ) 

I 

,\U4=N/4 
I.J4P1=NU4+1 
Nu4P2=ND4PJL  + l 
f U2P2=Hl<4+ND4P2 
UO  8 Lu=l»M 
LMX=2**(M-L0) 

LIX=2*L'-.X 

ISCL=N/lIX 

L j 8 LM=1»LMX 

lAKG=(LM-i)*IbCL+l 

IF (IARG.LE,N04P1)  GO  TO  4 

C--CC(ND2P2-IARG) 

b=IbM*CC(IARG-ND4) 

CO  TO  6 
4 C=CC(IARG) 

S=lSN*CC(ND4P2~IAKb) 
o i,0  8 LI=LIX#N»LIX 
J1=LI-LIX+LM 
J2=Jl+LMX 
T1=X( Jl)-X{ J2) 

T2=Y(Ul)-Y(J2) 

X(Jl)=A( J1)+X(J2) 

Y v vJl ) = Y ( J1)+Y(J2) 
x(J2)=C*T1-S*T2 

Y ^ J2)=C*T2+S*T 1 
ti  CONTINUE 

LU  40  Jzl » 12 
L(J)=1 

IF(J-M)  31#31»40 
31  UJ>=2**(M+1-J) 
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40  CONTINUE 
JN=1 

bO  60  J1=1#L1 
00  60  J2=J1»L2»LI 
bo  60  J3=vJ2,L3»L2 
CO  60  J4=U3rL4»L3 
bQ  60  J5=J4,Lb»L4 
DO  6C  J6-0brL6fLb 
bO  bO  J7=06#L7»L6 
bO  60  J8=J7»L6»L7 
bO  bO  09=J6»L9»L8 
LO  bO  J10-J9#L10#Ly 
DO  60  J11=J10»L11»U10 
00  60  JR=JH,L12»Lll 
Ih  ( JN-oR)  51,51,52 

51  K=a(JN) 

MJN)=X(JR) 

X i JR ) -K 
F1=YCJW) 

Y ( J|J)=Y  (JR) 

Y(OR)=FI 

52  wNSJN-fl 
60  CONTINUE 

KlTURN 

t-Nb 


SUbROUTINL  OTRCOS(CrN) 
bihcNSION  CU) 
N4i=N/4+l 
SCL=b.283185307/N 
bO  1 1=1, N41 

1 C(i)=COS( (I-1)*SCL) 
RETURN 
ENU 
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